Internal states of model isotropic granular packings. 
I. Assembling process, geometry and contact networks. 
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This is the first paper of a series of three, in which we report on numerical simulation studies 
of geometric and mechanical properties of static assemblies of spherical beads under an isotropic 
pressure. The influence of various assembling processes on packing microstructures is investigated. It 
is accurately checked that frictionless systems assemble in the unique random close packing (RCP) 
state in the low pressure limit if the compression process is fast enough, higher solid fractions 
corresponding to more ordered configurations with traces of crystallization. Specific properties 
directly related to isostaticity of the force-carrying structure in the rigid limit are discussed. With 
frictional grains, different preparation procedures result in quite different inner structures that 
cannot be classified by the sole density. If partly or completely lubricated they will assemble like 
frictionless ones, approaching the RCP solid fraction "I>rcp — 0.639 with a high coordination number; 
2* ~ 6 on the force-carrying backbone. If compressed with a realistic coefficient of friction fi — 0.3 
packings stabilize in a loose state with $ ~ 0.593 and z* ~ 4.5. And, more surprisingly, an idealized 
"vibration" procedure, which maintains an agitated, coUisional regime up to high densities results 
in equally small values oi z* while "1> is close to the maximum value ^rcp. Low coordination 
packings have a large proportion (>10%) of rattlers - grains carrying no force - the effect of which 
should be accounted for on studying position correlations, and also contain a small proportion of 
localized "floppy modes" associated with divalent grains. Low pressure states of frictional packings 
retain a flnite level of force indeterminacy even when assembled with the slowest compression rates 
simulated, except in the case when the friction coefficient tends to infinity. Different microstructures 
are characterized in terms of near neighbor correlations on various scales, and some comparisons with 
available laboratory data are reported, although values of contact coordination numbers apparently 
remain experimentally inaccessible. 

PACS numbers: 45.70.-n, 83.80.Fg, 46.65.+g, 62.20.Fe 



I. INTRODUCTION 
A. Context and motivations 

The mechanical properties of solidlike granular pack- 
ings and their microscopic, grain-level origins are an ac- 
tive field of research in material science and condensed- 
matter physics [J [3i 9 • Motivations are practical, origi- 
nated in soil mechanics and material processing, as well 
as theoretical, as general approaches to the rheology of 
different physical systems made of particle assemblies out 
of thermal equilibrium Q are attempted. 

The packing of equal-sized spherical balls is a simple 
model for which there is a long tradition of geometric 
characterization studies. Packings are usually classified 
according to their density or solid volume fraction $, and 
the frequency of occurrence of some local patterns. Di- 
rect observations of packing microstructure is difficult, 
although it has recently benefitted from powerful imag- 
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ing techniques [1, H, [3|- The concept of random close 
packing (RCP), is often invoked [ajsl, although some 
authors criticized it as ill-defined [10] ■ It corresponds 
to the common observation that bead packings without 
any trace of crystalline order do not exceed a maximum 
density, ^rcp, shghtly below 0.64 

Mechanical studies in the laboratory have been per- 
formed on granular materials for decades in the realm 
of soil mechanics, and the importance of packing frac- 
tion $ on the rheological behavior has long been rec- 
ognized [U, [nl, [H, The anisotropy of the packing 
microstructure, due to the assembling process, has also 
been investigated [H, [lB|, and shown to influence the 
stress-strain behavior of test samples [1^1 , as well as the 
stress field and the response to perturbations of gravity- 
stabilized sandpiles or granular layers [l7j . 

Discrete numerical simulation (isj proved a valuable 
tool to investigate the internal state of packings, as it 
is able to reproduce mechanical behaviors, and to iden- 
tify relevant variables other than <I>, such as coordination 
number and fabric (or distribution of contact orienta- 
tions) [H, [lO, mi, m, [2^ . In the case of sphere packings, 
simulations have been used to characterize the geometry 
of gravity-deposited systems [13, HH] or oedometrically 
compressed ones [2^ , to investigate the quasistatic, hys - 
teretic stress-strain dependence in solid packings [23, |23| , 
and their pressure-dependent elastic moduli in a com- 
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pression experiment [H, [s^l ■ 

However, in spite of recent progress, quite a few ba- 
sic questions remain unresolved. It is not obvious how 
closely the samples used in numerical simulations actu- 
ally resemble laboratory ones, for which density is often 
the only available state parameter. Both simulations and 
experiments resort to certain preparation procedures to 
assemble granular packings, which, although their influ- 
ence is recognized as important, are seldom studied, or 
even specified. One method to produce dense packings 
in simulations is to set the coefficient of intergranular 
friction to zero [1^, [l^l or to a low value 28;] in the 
assembling stage, while a granular gas gets compacted 
and equilibrated under pressure. On keeping friction- 
less contacts, this results, in the limit of low confin- 
ing stress, in dense systems with rather specific proper- 
ties [20, |31| , related to isostaticity and potential energy 
minimization [33 |. Examples of traditional procedures 
in soil mechanics are rain deposition under gravity, also 
known as air pluviation (which produces homo gen eous 
states if grain flow rate and height of free fall [33] are 
maintained constant) and layerwise deposition and dry 
or moist tamping. Those two methods were observed to 
produce, in the case of loose sands, different structures 
for the same packing density [U . Densely packed parti- 
cle assemblies can also be obtained in the laboratory by 
vibration, or application of repeated "taps" [H, [s^ to a 
loose deposit. How close are dense experimental sphere 
packings to model configurations obtained on simulating 
frictionless particles ? How do micromechanical parame- 
ters influence the packing structure ? Is the low pressure 
limit singular in laboratory grain packings and in what 
sense ? 



B. Outline of the present study 

The present paper provides some answers to such ques- 
tions, from numerical simulations in the simple case 
of isotropically assembled and compressed homogeneous 
packings of spherical particles. It is the first one of a 
series of three, and deals with the geometric characteri- 
zation of low pressure isotropic states assembled by dif- 
ferent procedures, both without and with intergranular 
friction. The other two, hereafter referred to as papers 
II 37] and III [11], respectively investigate the effects 
of compressions and pressure cycles, and the elastic re- 
sponse of the different numerical packings, with compar- 
isons to experimental results. Although mechanical as- 
pects are hardly dealt with in the present paper, we insist 
that geometry and mechanics are strongly and mutually 
related. We focus here on the variability of the coordi- 
nation numbers, which will prove important for mechan- 
ical response properties of granular packings, and show 
that equilibrated packs of identical beads can have a rel- 
atively large numbers of "rattler" grains, which do not 
participate in force transmission. We investigate the de- 
pendence of initial states on the assembling procedure. 



both with and without friction. We study the effects 
of procedures designed to produce dense states (close to 
RCP), and we characterize the geometry of such states 
on different scales. 

It should be emphasized that we do not claim here to 
mimic experimental assembling procedures very closely. 
Rather, we investigate the results of several prepara- 
tion methods, which are computationnally convenient, 
maintain isotropy, and produce equilibrated samples with 
rather different characteristics. Those methods neverthe- 
less share some important features with laboratory pro- 
cedures, and we shall argue that the resulting states are 
plausible models for experimental samples. 

The numerical model and the simulation procedures 
(geometric and mechanical parameters, contact law, 
boundary conditions) are presented in Section [TTl where 
some basic definitions and mechanical properties per- 
taining to granular packs are also presented or recalled. 
Part mil discusses the properties of frictionless packings, 
and introduces several characterization approaches used 
in the general case as well. Section HVl then describes dif- 
ferent assembling procedures of frictional packings and 
the resulting microstructures. Section |V| discusses per- 
spectives to the present study, some of which are pur- 
sued in papers II and HI of the series. Appendices deal 
with technical issues, and also present a more detailed 
comparison with some experimental data. 

This being a long paper, it might be helpful to specify 
which parts can be read independently. On first going 
through the paper, the reader might skip Section fill D[ 
dealing with a rather specific issue. The properties stated 
or recalled in Section Hi Cl are used to discuss stability is- 
sues and isostatic values of coordination numbers, but 
they can also be overlooked in a first approach. Finally, 
Section IIVI can be read independently from Section IIIII 
apart from the explanations about equilibrium conditions 
(in paragraph llll B 2| l and the treatment of rattlers (para- 
graph [TTrETl) ■ Sections IIIII and HVl both have conclusive 
subsections which summarize the essential results. 



II. MODEL, NUMERICAL PROCEDURES, 
BASIC DEFINITIONS 

A. Intergranular forces 

We consider spherical beads of diameter a (the value 
of which, as we ignore gravity, will prove irrelevant), in- 
teracting in their contacts by the Hertz law, relating the 
normal force N to the elastic normal defiection h as : 

N = ^h^/\ (1) 
3 ^ ' 

In Eqn. (TJ wc introduced the notation 
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E is the Young modulus of the beads, and u the Poisson 
ratio. For spheres, h, the elastic deflection of the contact, 
is simply the distance of approach of the centers beyond 
the first contact. The normal stiffness Kn of the contact 
is defined as the rate of change of the force with normal 
displacement: 

rfAT ^ :^/jl/2 ^ ^^2/3^1/3^1/3 (2) 

dh 2 2 ^ ' 

Although many geometric features of particle packings do 
not depend on the details of the model for contact elastic- 
ity, and could be observed as well with a simpler, linear 
unilateral elastic model, it is necessary to implement suit- 
able non-linear contact models to deal with the mechan- 
ical properties in papers II and III [s^, . Tangential 
elasticity and friction in contacts are appropriately de- 
scribed by the Cattaneo-Mindlin-Deresiewicz laws [soj . 
which we implement in a simplified form, as used e.g., in 
refs. [l^, im : the tangential stiffness Kt relating, in the 
elastic regime, the increment of tangential reaction dT to 
the relative tangential displacement increment dux is a 
function of h (or N) alone {i.e., it is kept constant, equal 
to its value for T = 0): 

dT — RxduT, with 

2-2iy ~ ^ , ,„ (3) 

Kt = Kn - EV^h^/^ 

2 — V 2 — V 

To enforce the Coulomb condition with friction coeffi- 
cient \x, T has to be projected back onto the circle of 
radius ^A'' in the tangential plane whenever the incre- 
ment given by eqn. [3] would cause its magnitude to ex- 
ceed this limit. Moreover, when A decreases to A — 8N , 
T is scaled down to the value it would have had if iV 
had constantly been equal to iV — in the past. It 
is not scaled up when N increases. Such a procedure, 
suggested e.g., in [i^l, avoids spurious increases of elas- 
tic energy for certain loading histories. More details are 
given in Appendix [XI 

Finally, tangential contact forces have to follow the ma- 
terial motion. Their magnitudes are assumed here not to 
be affected by rolling (i.e., rotation about a tangential 
axis) or pivoting (i.e., rotation about the normal axis), 
while their direction rotates with the normal vector due 
to rolling, and spins around it with the average spinning 
rate of the two spheres (to ensure objectivity). The cor- 
responding equations are given in Appendix [Bl 

In addition to the contact forces specified above, we 
introduce viscous ones, which oppose the normal rela- 
tive displacements (we use the convention that positive 
normal forces are repulsive): 

A^" = oi{h)h (4) 

The damping coefficient a depends on h, and we choose 
its value as a fixed fraction Q of the critical damping 
coefficient of the normal (linear) spring of stiffness Kf^(\i) 
(as given by ([2])) joining two beads of mass m: 

a(h) = CV^mK^^. (5) 



From (l2|), a is thus proportional to /i^/*, or to N^/^. 
The same damping law was used in [4l| . Admittedly, 
the dissipation given by (HI)-® has little physical justifi- 
cation, and is rather motivated by computational conve- 
nience. We shall therefore assess the influence of C on the 
numerical results. The present study being focussed on 
statics, we generally use a strong dissipation, C — 0.98, 
to approach equilibrium faster. This particular value is 
admittedly rather arbitrary: the initial motivation for 
choosing < 1 is the computational inefficiency of over- 
damped contacts with C > 1 in the case of linear contact 
elasticity. Yet we did not check whether values of 1 or 
even higher would cause any problem with Hertzian con- 
tacts. In the linear case, the restitution coefficient in a 
binary collision varies as a very lastly decreasing func- 
tion of Q, and changes of Q in the range between 0.7 and 
1 have virtually no detectable effect. 

We do not introduce any tangential viscous force, and 
impose the Coulomb inequality to elastic force compo- 
nents only. We choose the elastic parameters E = 70 GPa 
and ly = 0.3, suitable for glass beads, and the friction co- 
efficient is attributed a moderate, plausible value fj, = 0.3. 
These choices are motivated by comparisons to experi- 
mental measurements of elastic moduli, to be carried out 
in paper III [38l |. 

B. Boundary conditions and stress control 

The numerical results presented below were obtained 
on samples of n = 4000 beads, enclosed in a cubic or 
parallelipipedic cell with periodic boundary conditions. 

It is often in our opinion more convenient to use pres- 
sure (or stress) than density (or strain) as a control pa- 
rameter (a point we discuss below in Section IIIip . We 
therefore use a stress-controlled procedure in our sim- 
ulations, which is adapted from the Parrinello-Rahman 
molecular dynamics (MD) scheme [43 |. The simulation 
cell has a rectangular parallelipipedic shape with lengths 
La parallel to coordinate axes a (1 < a < 3). La values 
might vary, so that the system has 6A^-|-3 configurational 
degrees of freedom,which are the positions and orienta- 
tions of the A^ particles and lengths La. fl = L1L2L3 
denotes the sample volume. We seek equilibrium states 
with set values (Sa)i<Q;<3 of all three principal stresses 
o'aa ■ We use the convention that compressive stresses are 
positive. 

It is convenient to write position vectors r^, defining a 
square 3x3 matrix with L^'s on the diagonal, as 

{l<i<N) r, ^L-s„ 

Si denoting corresponding vectors in a cubic box of unit 
edge length. In addition to particle angular and linear 
velocities, which read 

{l<i<N) Vi = L-s» + t-Si, 

one should evaluate time derivatives L.^. Equations of 
motion are written for particles in the standard form. 
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i.e., (Fi denoting the total force exerted on grain i) 

{l<i<N) mii^=hr^ (6) 

and the usual equation for angular momentum. Mean- 
while, lengths La satisfy the following equation of mo- 
tion, in which Tjj is the vector joining the center of i to 
the center of j, subject to the usual nearest image con- 
vention of periodic boundary conditions: 



MLa = 



LIE' 



i<j 



(7) 



Ley 



Within square brackets on the right-hand side of Eqn. [7l 
one recognizes the familiar formula [l^ IH, Q for ilaaa , 
g_ being the average stress in the sample: 



C. Rigidity and stiffness matrices 

We introduce here the appropriate formalism and state 
the relevant properties of static contact networks. It is 
implied throughout this section that small displacements 
about an equilibrium configurations are dealt with to first 
order (as an infinitesimal motion, i.e. just like veloci- 
ties), and related to small increments of applied forces, 
moments and stresses. In the following we shall exploit 
the definitions of stiffness matrices K^^'' (Eqn. [T5)) and 

K^^^ to discuss stability properties of packings. The 
corrections to the degree of force indeterminacy due to 
free mechanism motions, as expressed by relations (jl9p 
or (PU)) . will also be used. 

The properties are stated in a suitable form to the peri- 
odic boundary conditions with controlled diagonal stress 
components, as used in our numerical study. 



i<j 



(8) 



All three diagonal stress components should thus equate 
the prescribed values at equilibrium. The acceleration 
term will cause the cell to expand in the corresponding 
direction if the stress is too high, and to shrink if it is 
too low. Eqn. [7] involves a generalized mass M associated 
with the changes of shape of the simulation cell. M is set 
to a value of the order of the total mass of all particles in 
the sample. This choice was observed to result in collec- 
tive degrees of freedom La approaching their equilibrium 
values under prescribed stress S somewhat more slowly 
(but not exceedingly so) than (rescaled) positions s. 

The original Parrinello-Rahman method was designed 
for conservative molecular systems, in such a way that 
the set of equations is cast in Lagrangian form. This 
implies in particular additional terms in ([6]), involving 
L. Such terms were observed to have a negligible influ- 
ence on our calculations and were consequently omitted. 
Granular materials are dissipative, and energy conserva- 
tion is not an issue (except for some elastic properties, see 
paper III 38]). Further discussion of the stress-controlled 
method is provided in Appendix [Cl 

Equations (O and ([7]), with global degrees of free- 
dom La slower than particle positions, lead to dynam- 
ics similar to those of a commonly used procedure in 
granular simulation (45j . This method consists in re- 
peatedly changing the dimensions of the cell by very 
small amounts, then computing the motion of the grains 
for some interval of time. A "servo mechanism" can be 
used to impose stresses rather than strains j30| . Our ap- 
proach might represent a simplification, as it avoids such 
a two-stage procedure. It should be kept in mind that 
we restricted our use of equation ([7]) to situations when 
changes in the dimensions of the simulation cell are very 
slow and gradual. The perturbation introduced in the 
motion of the grains, in comparison to the more familiar 
case of a fixed container, is very small. 



i. Defimtion of stiffness matrix 

We consider a given configuration with bead center 
positions (r^, 1 < i < n) and orientations {0i, 1 < i < n), 
and cell dimensions (La, 1 < a < 3). The grain center 
displacements (ui)i<i<„ are conveniently written as 



with a set of displacements satisfying periodic bound- 
ary conditions in the cell with the current dimensions, 
and the elements of the diagonal strain matrix e express 
the relative shrinking deformation along each direction, 
Ea = —ALa/ La. Gathering all coordinates of particle 
(periodic) displacements and rotation increments, along 
with strain parameters, one defines a displacement vector 
in a space with dimension equal to the number of degrees 
of freedom Nf — 6n + 3, 



U — ((ui, A6'i)i<i<„, (eQ)i<a<3) . 



(9) 



Let A'';; denote the number of intergranular contacts. In 
every contacting pair i-j, we arbitrarily choose a "first" 
grain i and a "second" one j. The normal unit vector 
iiij points from i to j (along the line joining centers for 
spheres). The relative displacement Suij is defined for 
spherical grains with radius R as 



X Ruij — iij 



9j X Ruij 



(10) 



in which is the vector pointing from the center of 
the first sphere i to the nearest image of the center of 
the second one j. The normal part Sufj of Suij is the 
increment of normal deflection hij in the contact. (jlOp 
defines a 3A''c x Nf matrix G which transforms U into 
the 3Ac-dimensional vector of relative displacements at 
contacts Su: 



Su^G V 



(11) 
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In agreement with the hterature on rigidity theory of 
frameworks [i^ {—G_ is termed normahzed rigidity ma- 
trix in that reference) , we call G the rigidity matrix. 

In each contact a force Fij is transmitted from i to j, 
which is split into its normal and tangential components 
as Fij = NijUij + T.y . The static contact law (without 
viscous terms) expressed in Eqns. ([IJ, Q, with the condi- 
tions stated in Section Hi Al relates the 3A^c-dimensional 
contact force increment vector Af , formed with the val- 
ues AiVy , ATy- of the normal and tangential parts of all 
contact force increments, to 6u: 



At^lC- 5u. 



This defines the {3Nc x 3Nc) matrix of contact stiffnesses 
IC. JC is block diagonal (it does not couple different con- 
tacts), and is conveniently written on using coordinates 
with YYij as the first basis unit vector. In simple cases the 
3x3 block of ^ corresponding to contact i, j, JC.. is di- 
agonal itself and contains stiffnesses Kpf(hij) and (twice 
in 3 dimensions) Kxihij) as given by ^ and 



fC.. ^ 



KN{h,j) 
Kt{K,) 
KT{h,j) 



(13) 



More complicated non-diagonal forms of 1C_^. , which actu- 
ally depend on the direction of the increments of relatives 
displacements in the contact, arc found if friction is fully 
mobilized (which does not happen in well-equilibrated 
configurations) , or corresponding to those small motions 
reducing the normal contact force. The effects of such 
terms is small, with our choice of parameters, and is dis- 
cussed in paper III (38] . 

External forces F^ and moments Ti (at the center) ap- 
plied to the grains, and diagonal Cauchy stress compo- 
nents Sq, can be gathered in one A^j^-dimensional load 
vector F™*: 



((F„r,) {VlY. a) l<a<:i) , 



(14) 



chosen such that the work in a small motion is equal to 
pcxt The equilibrium equations - the statements that 
contact forces f balance load F"'^* - is simply written with 
the tranposed rigidity matrix, as 



= . f . 



(15) 



This is of course easily checked on writing down all force 
and moment coordinates, as well as the equilibrium form 
of stresses: 



(16) 



1<J 



As an example, matrices G and ^G were written down 
in [13] in the simple case of one mobile disk with 2 con- 
tacts with fixed objects in 2 dimensions, the authors re- 
ferring to — ^G as the "contact matrix" . The same defi- 
nitions and matrices are used in ^A§\ in the more general 
case of a packing of disks. 



Returning to the case of small displacements associated 
with a load increment AF°^*, one may write, to first 
order in U, 



(17) 



with a total stiffness matrix K, comprising two parts, 

K^^-* and K^^-*, which we shall respectively refer to as 

the constitutive and geometric stiffness matrices. K 
results from Eqns. [TTl (T^l and [T5l 

K^^' = '^G-IC-G. (18) 



(12) K'^' is due to the change of the geometry of the pack- 



ing. Its elements (see Appendix [B|) . relatively to to their 
counterparts in K*-^-*, are of order F/KnR ^ h/R, and 
therefore considerably smaller in all practical cases. The 
constitutive stiffness matrix is also called "dynamical ma- 
trix" [m im . One advantage of decomposition (fT5|) is to 
separate out the effects of the contact constitutive law, 
contained in IC and those of the contact network, con- 
tained in G. G is sensitive in general to the orientations 
of normal unit vectors Uij and to the "branch vectors" 
joining the grain centers to contact positions - which re- 
duce to Rriij for spheres of radius R. K*-^\ on the other 
hand, unlike G, is sensitive to the curvature of grain sur- 
faces at the contact point [i^, [H^] • 



2. Properties of the rigidity matrix 

To the rigidity matrix are associated the concepts 
(familiar in structural mechanics) of force and velocity 
(or displacement) indeterminacy, of relative displacement 
compatibility and of static admissibility of contact forces. 
Definitions and properties stated in [32] for frictionless 
grains, straightforwardly generalize to packings with fric- 
tion. 

The degree of displacement indeterminacy (also called 
degree of hypostaticity [13]) is the dimension k of the 
kernel of G, the elements of which are displacements vec- 
tors U which do not create relative displacements in the 
contacts: Su = 0. Such displacements are termed (first- 
order) mechanisms. Depending on boundary conditions, 
a grain packing might have a small number kg of "triv- 
ial" mechanisms, for which the whole system moves as 
one rigid body. In our case, attributing common values 
of u to all grains gives ko = 3 independent global rigid 
motions. 

The degree of force indeterminacy H (also called degree 
of hyperstaticity [s^]) is the dimension of the kernel of 
^G. or the number of independent self-balanced contact 
force vectors. If the coordinates of f are regarded as 
the unknowns in system of equations (fTS]) . and if F'"^* 
is supportable, then there exists a whole H-dimensional 
affine space of solutions. 

From elementary theorems in linear algebra one de- 
duces a general relation between H and k |32| 



Nf + -ii^3Nc + k. 



(19) 
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An isostatic packing is defined as one devoid of force 
and velocity indeterminacy (apart from trivial mecha- 
nisms). Excluding trivial mechanisms (thus reducing Nf 
to Nf — ko), and loads that are not orthogonal to them, 
one then has a square, invertible rigidity matrix. To 
any load corresponds a unique set of equilibrium contact 
forces. To any vector of relative contact displacements 
corresponds a unique displacement vector. 

With frictionless objects, in which contacts only carry 
normal forces, it is appropriate to use A^c-dimensional 
contact force and relative displacement vectors, contain- 
ing only normal components, and to define the rigidity 
matrix accordingly [321 . Then should be written as 



whence an exact relation between P and contact deflec- 
tions: 



Nf Nc 



(20) 



In the case of frictionless spherical particles, all rotations 
are mechanisms, hence a contribution of 3 to k. Thus 
one may in addition ignore all rotations, and subtract 
3n both from Nf and from k, so that ([^0)1 is still valid. 
In such a case, the rigidity matrix coincides (up to a 
sign convention and normalization of its elements) with 
the one introduced in central-force networks, trusses and 
tensegrity structures [5l[. Donev et al, in a recent pub- 
lication on sphere packings f53| , call rigidity matrix what 
we defined as its transpose "^G. 



D. Control parameters 

The geometry and the mechanical properties of sphere 
packings under given pressure P depends on a small set 
of control parameters, which can be conveniently defined 
in dimensionless form [23, ^] . 

Such parameters include friction coefficient /i and vis- 
cous dissipation parameter Q, which were introduced in 
SecHEl 

The elastic contact law introduces a dimensionless 
stiffness parameter k, which we define as: 




(21) 



Note that n does not depend on bead diameter a. Under 
pressure P, the typical force in a contact is of order Pa^ . 
It corresponds to a normal deflection h such that Pa^ ~ 
E^/ah^/"^ due to the Hertz law ([T]). Therefore, k sets 
the scale of the typical normal deflection h in Hertzian 
contacts, as h/a ^ 1/k. 

In the case of monodisperse sphere packings in equi- 
librium in uniform state of stress a, pressure P = ira_/i 
is directly related to the average normal force {N). Let 
us denote as $ the solid fraction and z the coordination 
number (z = 2Nc/n). As a simple consequence of the 
classical formula for stresses recalled in Sec. Ill Bl (Eqn. [5] 
in the static case, or Eqn. fTH)) . one has, neglecting contact 
deflections before diameter a, 



(fe3/2) 
a3/2 



The limit of rigid grains is approached as k ^ oo. k 
can reach very high values for samples under their own 
weight, but most laboratory results correspond to levels 
of confining pressure in the 100 kPa range. Experimental 
data on the mechanical properties of granular materials 
in quasistatic conditions below a few tens of kPa are very 
scarce (see, however, 0] and [l^). This is motivated by 
engineering applications (100 kPa is the pressure below a 
few meters of earth) , and this also results from difficulties 
with low confining stresses. Below this pressure range, 
stress fields are no longer uniform, due to the influence of 
the sample weight, and measurements are difficult {e.g., 
elastic waves of measurable amplitude are very strongly 
damped) . 

We set the lowest pressure level for our simulation of 
glass beads to 1 kPa or 10 kPa, which corresponds to 
K ~ 181000 and k ~ 39000. Such values, as we shall 
check, are high enough for some characteristic proper- 
ties of rigid sphere packings to be approached with good 
accuracy. Upon increasing P, the entire experimental 
pressure range will be explored in the two companion 
papers [13, [11]. 

Another parameter associated with contact elasticity 
is the ratio of tangential to normal stiffnesses (constant 
in our model), related to the Poisson ratio of the material 
the grains are made of. Although we did not investi gate 
the role of this parameter, several numerical studies |4ll . 
[56| showed its influence on global properties to be very 
small. 

The "mass" M of the global degrees of freedom is cho- 
sen to ensure slow and gradual changes in cell dimensions, 
and dynamical effects are consequently assessed on com- 
paring the strain rate e to intrinsic inertial times, such as 
the time needed for a particle of mass m, initially at rest, 
accelerated by a typical force Pa^, to move on a distance 
a. This leads to the deflnition of a dimensionless inertia 
parameter : 



I = t\pmfaP . 



(23) 



P 



(22) 



The quasistatic limit can be defined as / — s- 0. / was 
successfully used as a control parameter in dense granu- 
lar shear flows [s^, [sl, [5§| , which might be modelled on 
writing down the / dependence of internal friction and 
density [63,[6l|. 

The sensitivity to dynamical parameters / and C 
should be larger in the assembling stage (as studied in the 
present paper) than in the subsequent isotropic compres- 
sion of solid samples studied in paper II [37(], for which 
one attempts to approach the quasi-static limit. In the 
following we will assess the influence of parameters ^, k, 
/ and C, on sample states and properties. 
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III. LOW-PRESSURE ISOTROPIC STATES OF 
FRICTIONLESS PACKINGS. 

A. Motivations 

Numerical samples are most often produced by com- 
pression of an initially loose configuration (a granular 
gas) in which the grains do not touch. If the friction 
coefficient is set to zero at this stage, one obtains dense 
samples, which depend very little on chosen mechani- 
cal parameters. These frictionless configurations are in a 
particular reference state which was recently investigated 
by several groups [3l|, [s^] . We shall dwell on such an aca- 
demic model as assemblies of rigid or slightly deformable 
frictionless spheres in mechanical equilibrium for several 
reasons. First, we have to introduce various characteri- 
zations of the microstructure of sphere packings that will 
be useful in the presence of friction too. Then, such sys- 
tems possess rather specific properties, which are worth 
recalling in order to assess whether some of them could 
be of relevance in the general case. Frictionless packings 
also represent, as we shall explain, an interesting limit 
case. Finally, one of our objectives is to establish the 
basic uniqueness, in the statistical sense, of the internal 
state of such packings under isotropic, uniform pressures, 
provided crystallization is thwarted by a fast enough dis- 
sipation of kinetic energy. 

B. Assembling procedures 

1. Previous results 

Since we wish to discuss a uniqueness property, we 
shall compare our results to published ones whenever 
they are available. Specifically, we shall repeatedly refer 
to the works of O'Hern, Silbert, Liu and Nagel [3l|, and 
of Donev, Torquato and Stillinger [52*1, hereafter respec- 
tively abbreviated as OSLN and DTS. Both are numeri- 
cal studies of frictionless sphere packings under isotropic 
pressures. 

OSLN use elastic spheres, with either Hertzian or lin- 
ear contact elasticity. They control the solid fraction $, 
and record the pressure at equilibrium. Their samples 
(from a few tens to about 1000 spheres) are requested to 
minimize elastic energy at constant density. For each one, 
pressure and elastic energy vanish below a certain thresh- 
old packing fraction <I>o, which is identified to the classical 
random close packing density. Above <f>o, pressure and 
elastic constants are growing functions of density. OSLN 
report several power law dependences of geometric and 
mechanical properties on $ — $0 which we shall partly 
review. 

DTS differ in their approach, as unlike OSLN (and 
unlike us) they use strictly rigid spherical balls, and 
approach the density of equilibrated rigid, frictionless 
sphere packing from below. They use a variant of the 
classical (event-driven) hard-sphere molecular dynamics 



method |44l . |62| | , in which sphere diameters are contin- 
uously growing, the Lubachevsky-Stillinger (LS) algo- 
rithm [sl, [ill, to compress the samples. DTS's approxi- 
mation of the strictly rigid sphere packing as the limit of 
a hard sphere glass with very narrow interstices (gaps) 
between colliding neighbors (contact forces in the static 
packing are then replaced by transfers of momentum be- 
tween neighbors), and their resorting to linear optimiza- 
tion methods [32 . [65j , enable them to obtain very accu- 
rate results in samples of 1000 and 10000 beads. 

DTS expressed doubts as to whether numerical "soft" 
(elastic) sphere systems could approach the ideal rigid 
packing properties, and both groups differ in their actual 
definition of jamming and on the relevance and definition 
of the random close packing concept. Relying on our own 
simulation results, we shall briefly discuss those issues in 
the following. 



2. Frictionless samples obtained by MD 

Our numerical results on packings assembled with- 
out friction are based on five different configurations of 
n — 4000 beads prepared by compression of a granular 
gas without friction. First, spheres are placed on the 
sites of an FCC lattice at packing fraction $ = 0.45 (be- 
low the freezing density, $ ~ 0.49 [66]). Then they are 
set in motion with random velocities, and left to inter- 
act in collisions that preserve kinetic energy, just like 
the molecules of the hard-sphere model fluid studied in 
liquid state theory [H, [11] . We use the traditional event- 
driven method [63], in a cubic cell of fixed size, until the 
initial crystalline arrangement has melted. Then, veloci- 
ties are set to zero, and the molecular dynamics method 
of Section |TT] is implemented with an external pressure 
equal to 10 kPa for glass beads (k ~ 39000). Energy is 
dissipated thanks to viscous forces in contacts, and the 
packing approaches an equilibrium state. Calculations 
are stopped when the net elastic force on each particle 
is below lO^^'a^P, the elastic contributions to the stress 
components equal the prescribed value P with relative er- 
ror smaller than 10~^ and the kinetic energy per particle 
is below lO^^Pa^. On setting all velocities to zero, it is 
observed that the sample does not regain kinetic energy 
beyond that value, while the unbalanced force level does 
not increase. We have thus a stable equilibrium state. 
This is further confirmed by the absence of mechanism in 
the force-carrying contact network, apart from the trivial 
free translational motion of the whole set of grains as one 
rigid body. From (jlSp , mechanisms coincide with "floppy 
modes" of the constitutive stiffness matrix K*-^-* (i.e., the 

elements of its kernel). The geometric stiffness K*-^\ as 
checked in Appendix [Bl is a very small correction (com- 
pared to those of K*-^-* the elements of matrix K^^-* are of 
order k~^). 

In the following, such configurations assembled with- 
out friction will be referred to as A states. 
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In order to check for a possible influence of the assem- 
bling procedure on the final configurations, we simulated 
another, similar sample series, denoted as A', for which 
the LS algorithm was used to bring the solid fraction from 
0.45 to 0.61, before equilibrating at the desired pressure 
with Hertzian sphere molecular dynamics. 

Observed geometric and mechanical characteristics of 
A and A' states are reported below and compared to 
other published results, in particuler those of OSLN and 
DTS. We also state specific properties of rigid frictionless 
sphere packings, to which A configurations at high n are 
close. Unlike OSLN, we use pressure or stiffness level k 
as the control parameter. The state OSLN refer to as 
"point J" , which appears as a rigidity threshold $ = $o 
if solid fraction $ is used as the control parameter, is 
approached here as /c — > c». 



3. Compression rates and duration of agitation stage 

Molecular dynamics is not the fastest conceivable route 
to minimize the sum of elastic and potential energies, and 
the MD approach does not necessarily find the nearest 
minimum in configuration space. For that purpose, the 
direct conjugate gradient minimization approach, as used 
by OSLN, which involves no inertia and follows a path 
of strictly decreasing energy in configuration space is the 
best candidate. 

However, the time scales involved in the MD simula- 
tions can be compared to experimental ones. In sim- 
ulations, A configurations approach their final density 
within a few tens of time r — y^m/ (aP), and come to 
their final equilibrium with a few hundreds of r. Com- 
parable laboratory experiments in which dense samples 
are assembled are sample preparation with the pluvia- 
tion or rain deposition technique, in which grains arc 
deposited at constant flow rate under gravity, with a con- 
stant height of free fall [H, [H, 113, HI] . Such an assem- 
bling technique produces homogeneous samples. Grains 
are first agitated near the free surface, and then subjected 
to a quasistatic pressure increase as pouring procedes. 
The relevant pressure scale corresponds therefore to the 
weight of the agitated superficial layer of the sample be- 
ing assembled [67|], typically of the order of 10 diameters, 
hence P ^ lOmg/a^ and t ~ y/a/{lOg), about 3 x 10^^ s 
for a = 1 mm. Approximating the compaction time by 
the time needed to renew entirely the agitated superficial 
layer, we obtain a few times 10~^ s if this time is to be 
of order lOr, as in our simulations. This corresponds to 
a fraction of a second to fill up a 10 cm high container, a 
value within the experimental range. The main conclu- 
sion from this crude analysis is that laboratory assem- 
bling processes are rather fast, with typical compaction 
times similar to those of our simulated isotropic compres- 
sion procedure. 

On the other hand, the LS procedure followed by DTS, 
which we used to produce our A' samples, unavoidably 
involves many collisions and a considerable level of agi- 



tation while particle diameters grow at a prescribed rate. 
In practice, kinetic energy actually increases on imple- 
menting the LS algorithm: receding velocities after a col- 
lision have to be artificially increased in order to make 
sure particles that are continuously growing in diame- 
ter actually move apart after colliding [g^ [g^] ■ Veloci- 
ties have to be scaled down now and then for computa- 
tional convenience, a feature the actual compacting pro- 
cess, depending on the ratio of growing rate to quadratic 
velocity average, is sensitive to. In our implementation 
there were typically 110 collisions per sphere in the range 
0.49 < $ < 0.58 (the most dangerous interval for crystal 
nucleation [fii,!!!]), and 90 collisions for 0.58 < $ < 0.61. 
DTS report using expansion rates of 10^^, while ours 
started as 10~^, in units of the quadratic mean velocity. 

Consequently, the order of simulation results, from the 
fastest, least agitated case to the slowest one is as follows: 
first the OSLN results, then our A, followed by our A' 
series, and finally the simulations by DTS (who used a 
slower implementation of the LS method than our A' 
one). 

C. Energy minimization and density 

1. What is "jamming" ? 

In spite of a long tradition of studies on the geometry 
of sphere assemblies, the connection between mechanical 
equilibrium and density maximization has seldom been 
stressed. This property was presented, in slightly differ- 
ent forms, in the mathematics [70] and physics [s^ . [s^ 
literature. It is worth recalling it here, as the purpose 
of this work is to discuss both geometric and mechan- 
ical properties of such particle packings. This connec- 
tion is simply expressed on noting that configurations of 
rigid, frictionless, non-adhesive spherical particles in sta- 
ble equilibrium under an isotropic confining pressure are 
those that realize a local minimum of volume in configu- 
ration space, under the constraint of mutual impenetra- 
bility. It is no wonder then that the isotropic compaction 
of frictionless balls is often used as a route to obtain 
dense samples [13, IsS] • In DTS [13] and in other works 
by the same group |65l . Fn| , the authors use a definition of 
strictly jammed configurations of hard particles as those 
for which particles cannot move without interpenetrating 
or increasing the volume of the whole system. Their def- 
inition is therefore exactly equivalent to that of a stable 
equilibrium state with rigid, frictionless grains under an 
isotropic confining pressure. 

If we now turn to elastic, rather than rigid, spherical 
particles, with Hertzian contacts as defined in Sec. m 
then stable equilibrium states under given pressure P 
are local minima of the potential energy defined as {H 
denotes the Heavside step function) 

W = Pn+ ^/.fiJ(/.,). (24) 

l<i<j <n 
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As stiffness parameter k increases, the second term 
of ([M]) imparts an increasing energetic cost to elastic 
deflections hij , and the solution becomes an approxima- 
tion to a minimum of the first term, with impenetrability 
constraints, i.e., a stable equilibrium state of rigid, fric- 
tionless balls. The value of k is an indicator of the dis- 
tance to the ideal, rigid particle configuration, and it is 
arguably more convenient to use that the density, used by 
OSLN, because it does not vary between different sam- 
ples. OSLN had to adjust the density separately for each 
sample in order to approach the limit of rigid grains, so 
that the pressure approached zero, corresponding to a 
rigidity threshold. Their definition of jamming is based 
on a local minimum of elastic energy, and therefore also 
coincides with ours: a jammed state is a stable equilib- 
rium state. 



2. Solid fractions 

Our A configurations have a solid fraction <!> = 
0.6370±0.0015 (indicated error bars correspond through- 
out the paper to one sample-to-sample standard devia- 
tion). We shall check below that the small density dif- 
ference between n = 39000 and k ^ cx3 is much smaller 
than the statistical uncertainty on <i>. OSLN performed 
a careful statistical analysis of finite size effects and un- 
certainty on leading to estimates shown on Fig. [T] 
Fig. [1] also shows another MD data point we obtained for 
n = 1372. Our $ values coincide with OSLN's estima- 
tion of size-dependent averages and fluctuations, once it 
is extrapolated to larger sample sizes (or, possibly, our 
configurations are very slightly denser). Our A' samples 
exhibit higher densities than A ones, $ = 0.6422 ±0.0002 
- a fairly small difference, but clearly larger than error 
bars. DTS do not report $ values very precisely, but 
mention solid fractions in the range 0.625 to 0.63 jH, 
page 7], on excluding the volume of rattlers, particles 
that transmit no force. This entails 0.639 < $ < 0.644 
once those inactive grains, which represent about 2.2% 
of the total number, are taken into account. LS-made 
samples were shown in [lol | to jam, depending on the 
compression rate, over the whole solid fraction range be- 
tween 0.64 and the maximum value 7r/(3\/2) correspond- 
ing to the perfect FCC crystal. The final values of the 
solid fraction therefore correlate with the duration of the 
agitation stage in the assembling procedure. The RCP 
density is traditionnally associated with a minimization 
of crystalline order. In the next section we check for indi- 
cations of incipient crystalline order in A and A' samples. 



D. Traces of crystalline order 

The possible presence of crystal nuclei, the FCC and 
HCP lattices (the former the more stable thermodynam- 
ically) and hybrids thereof being the densest possible ar- 
rangements, is a recurring issue in sphere packing studies. 
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FIG. 1: (Color online) Solid fraction $ versus n~ ' . Blue 
dotted line: average value, with standard deviation indicated 
as error bars, according to OSLN's results, extrapolated to 
71 ^ oo. Black round dots with error bars: our A samples for 
n = 4000 and other, very similar results for n = 1372. Square 
dot: A' samples with n — 4000 (this point has a smaller error 
bar). 



A recent numerical study of crystallization dynamics in 
the hard sphere fluid is that of Volkov et al. [6^ , in which 
the authors used several indicators and measures of in- 
cipient crystalline order that we apply here to A and A' 
states. First, bonds are defined as (fictitious) junctions 
between the centers of neighboring spheres if their dis- 
tance is smaller than some threshold, often chosen as 
corresponding to the first minimum in pair corelation 
function g{r) (about 1.4a in our case, see Sec. IIII F II) . 
Then, a local order parameter is associated to each grain 
i, as: 



47r 



rn—l 



n 1/2 



21 



m—~l 



(25) 



in which qim{i) is an average over all neighbors j of i 
numbered from 1 to Nh(i), the number of bonds of i: 



(26) 



nij denoting as usual the unit vector pointing from the 
center of i to the center of j. 

Qi and Qe, in particular, have been used to distin- 
guish different local orders 0, [sl, [6^ . In the sequel we 
use the average Qg"^ of (p6|) over all grains, as well as a 
global parameter Qg, defined on taking the average over 
all bonds within the sample, instead of those of a partic- 
ular grain i in (PS)) . The values of those parameters are 
given in table [D Global Qq values are small in large sam- 
ples, because they tend to average to zero in the presence 
of randomly oriented polycrystalline textures. They can 
be used nevertheless to observe crystallization in samples 
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of ~ 10000 beads, as they finally reach values comparable 
to the perfect crystal one [69| . 

Next, following [bB], we normalize the set 
{qim{i))-i<m<h on multiplying, for any given I and 
j, each of its 21 + 1 components by an appropriate 
common factor thus obtaining (jlimi^)) -i<in<u such that 



to the pair of values Q 



loc/ 



Q^°^{i)- We compared the 



m—l 

E 

m— — l 



(27) 



If the values qim{i) are viewed as the components of a 2Z + 
1-dimensional local order parameter, then qimii) niight 
be viewed as its "phase" or "angular" part, characteristic 
of the choice of a direction, rather than of the intensity 
or extent with which the system is locally ordered. Then 
a bond is termed crystalline if it joins two particles for 
which those "phases" are sufficiently correlated: (the star 
indicates complex conjugation) 



> 0.5. 



(28) 



A particle is said to be in a crystalline configuration if at 
least 7 of its bonds (out of 12.5-13, see table[l|) is "crys- 
talline" , according to definition (pS)) with 1 — 6. One may 
check how numerous those particles are and whether they 
tend to cluster in crystalline regions. Table |T] contains 
those various indicators, as observed in samples of type 
A and A' at the largest studied stiffness level. Order pa- 
rameters have a very small value, indicating as expected 
a large distance to crystal order. Only a small fraction 
of bonds and grains are declared "crystalline" according 
to the above definitions. However, it does transpire from 
the data of table U that A' states are consistently more 
"ordered" than A ones, with a small, but systematic dif- 
ference for all listed indicators (see also Appendix iPj) . 
Most notable is the increase of the size of "crystalline" 
regions. A direct visualization of those domains, as we 
checked, shows that they arc quite far from perfectly or- 
dered, but reveals some local tendency to organization 
in parallel stacked layers, and to the formation of 2D 
triangular lattice patterns within the layers. Luchnikov 
et al. [6§| report simulation of 16000 particle samples of 
the hard sphere fluid evolving towards crystallization at 
constant density (between 4> = 0.55 and $ — 0.6), as 
monitored by the global Qg parameter and the distribu- 
tion of local Qe values. They observed, like Volkov et 
al. [6^, that several thousands of collisions per particle 
were necessary for a significant evolution to take place, 
which is compatible with our observation of a detectable, 
but very small tendency with about 100 collisions per 
particle with our A' samples. 

Q^q'^ and 04°^^, as defined in (PS)) . were also used by 
Aste et al. [7| to characterize local arrangements, in an 
experimental study of sphere packing geometry by X-ray 
tomography. These results rely on observations of large 
samples of tens to hundreds of thousands of beads, al- 
though not isotropic. Particles are classified according 



geometry of our numerical samples of similar density to 
those experimental data, with the result that although 
the most frequently observed values of (5g°'^(i) and Q^^'^ii) 
were quite close to experimental ones in dense samples, 
and the proportion of hep- like particles were similar, fcc- 
like local environments were exceptional in simulations, 
whereas a few percent of the spheres were classified in 
that category in the experimental results. Quantitative 
results are given in Appendix [Dl Nucleation of crys- 
talline order is stron gly sensitive to sample history and 
boundary conditions [66| . 



E. Properties of force networks 

1. Identification and treatment of rattlers 

The rattlers are defined as the grains that do not 
participate in carrying forces and remain, therefore, 
free to "rattle" within the cage formed by their force- 
carrying, rigidly fixed neighbors. We refer to the net- 
work of contacting grains that carry forces as the back- 
bone. The backbone is the structure formed by non- 
rattler grains. The fraction xq of rattlers a,t k — 39000 is 
xo — 0.013± 0.002 in A samples, and it is slightly higher, 
xo = 0.018 ± 0.002 in A' ones. DTS report xq ^ 0.022, 
and hence once again our A' results are closer to theirs. 
The proportion of rattlers increases slightly for stiffer 
contacts (higher k values). 

Distinguishing between the backbone and the rattlers 
requires some care, as very small forces on the backbone 
might be confused with forces below tolerance between 
rattler and backbone grains, or between two rattlers. We 
apply the following simple procedure. First, we regard 
as rattlers all spheres having less than four contacts: less 
than three contacts implies a mechanism, and only three 
is impossible if forces are all strictly compressive. We also 
discard from the backbone all spheres with only forces 
smaller than the tolerance. Then, all the contacts of elim- 
inated spheres being also removed, other spheres might 
(although this is an extremely rare occurrence) have less 
than four contacts, so the procedure is iterated (twice at 
most is enough in our samples, although one such sweep 
is usually enough) until no more rattler is detected. We 
found this method to work correctly for n — 4000 and 
K — 39000. If one eliminates too many particles, the 
identified backbone might become floppy (hypostatic). 
We check, however, that its constitutive stiffness matrix 
remains positive definite, thereby avoiding such pitfall. 
The proportion of rattlers is likely to increase for stiffer 
contacts (higher k). 

The presence of rattlers complicates the analysis of ge- 
ometric properties of static packings, because their posi- 
tions are not determined by the equilibrium requirement. 
The rattlers are free to move within a "cage" formed by 
their backbone neighbors, and there is no obvious way in 
principle to prefer one or another of their infinitely many 
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State del a 


Z Zcr Qa 








QT 


^loCjCr 


Xcr 


(ricr) 


A 1.35 


12.36 ±0.03 2.91 ±0.06 (1.7 ±0.3) 


X 


10" 


-2 


0.392 ±0.001 


0.417 ±0.003 


0.080 ± 0.005 


19.8 


A' 1.35 


12.50 ±0.02 3.13 ±0.11 (1.9 ±0.5) 


X 


10' 


-2 


0.398 ± 0.0005 


0.420 ±0.002 


0.104 ±0.006 


54.8 


A 1.40 


13.11 ±0.02 2.94 ±0.06 (1.6 ±0.3) 


X 


10' 


-2 


0.370 ± 0.001 


0.394 ± 0.002 


0.083 ± 0.006 


22.8 


A' 1.40 


13.20 ±0.02 3.16 ±0.11 (1.8 ±0.3) 


X 


10' 


-2 


0.377 ± 0.0006 


0.397 ±0.003 


0.103 ±0.006 


64.5 



TABLE I: Indicators of possible incipient crystalline order in states A and A' at k = 39000. Z is the coordination number of 
first heighbors, Zcr the "crystalline bond" coordination number, Qg and Q^^"^ the global and (average) local order parameters, 
Qk>c,cr average value within crystalline regions, x^r the fraction of "crystalline" particles and (ricr) the mass average of the 
number of particles in a "crystalline cluster". First neighbors are defined here as those closer than distance dc = 1.35a or 
dc = 1.40a, near the first minimum in g{r). 



possible positions. This renders the evaluation of geo- 
metric data like pair correlations somewhat ambiguous. 
Moreover, rattlers, although scarce in frictionless pack- 
ings, can be considerably more numerous in frictional 
ones (see Section ITVT) . We therefore specify whether the 
results correspond to direct measurements on the con- 
figurations resuhing from the simulations, with rattlers 
floating in some positions resulting from compaction dy- 
namics, or whether rattlers have been fixed, each one 
having three contacts with the backbone (or some previ- 
ously fixed other rattler). To compute such fixed rattler 
positions with MD, we regard each backbone grain as 
a fixed object, exert small isotropically distributed ran- 
dom forces on all rattlers and let them move to a final 
equilibrium position (assuming frictionless contacts). A 
third possibility is to eliminate rattlers altogether before 
recording geometric data. These are three choices re- 
ferred to as I, II and III in the sequel, and we denote 
observed quantities with superscripts I, II or III accord- 
ingly 

Packings under gravity, if locally in an isotropic state 
of stress, are expected to be in the same internal state and 
to exhibit the same properties as the ones that are simu- 
lated here. In such a situation, individual grain weights 
are locally, within an approximately homogeneous sub- 
system, dominated by the externally imposed isotropic 
pressure. There is no rattler under gravity, but some 
grains are simply feeling their own weight, or perhaps 
that of one or a few other grains relying on them. Such 
grains are those that would be rattlers in the absence of 
gravity. Instead of freely floating within the cage of their 
backbone neighbors, they are supported by the cage floor. 
The situation should therefore be similar to that of our 
samples after all rattlers have been put in contact with 
the force-carrying structure (treatment II), except that 
the small external forces applied to them are all directed 
downwards. 



2. Coordination numbers 

Table [III gives the distribution of local coordination 
number values among the spheres for A and A' states 
at K = 39000. In this table, Xi is the proportion of 
grains with i contacts. If rattlers are stuck to the back- 



State 


a^o Xl X2 X3 X4, X5 Xq X7 Xs Xg Xio Xu 


A (I) 


1.3 11.1 23.2 28.4 22.6 10.3 2.8 0.3 0.02 


A' (I) 


1.8 11.6 22.5 28.2 22.3 10.6 2.8 0.3 0.01 


A (II) 


1.1 10.7 22.7 28.2 22.9 10.9 3.1 0.3 0.02 


A' (II) 


1.7 10.9 21.9 28.2 22.5 11.4 3.1 0.3 0.01 



TABLE 11: Percentages Xi of grains having i contacts in A and 
A' configurations, on ignoring the rare contacts with or be- 
tween rattlers (I), or on fixing the rattlers onto the backbone 
with small (randomly oriented) forces (II). 



bone (method II), one records slightly changed propor- 
tions of spheres with n > 3 contacts, to which values 
observed within samples under gravity should be com- 
pared. Distributions of local coordination numbers ob- 
served by DTS coincide to the data of Table Hi] within 
1%. We attribute this small difference to the influence 
of contact deflections of order K~^a in the MD results, 
while the DTS results are closer to ideally rigid packings 
(approached as open gaps tend to zero). 



3. Isostaticity 

We now discuss how the isostaticity property [3^ . [t^ . 
S III III, III, H III of equilibrium states of rigid, fric- 
tionless spheres, influences high k conflgurations of type 
A. 

Isostaticity is a property of the backbone, i.e. the 
force-carrying contact network, in equilibrium packings 
of rigid, frictionless spheres. It means that such networks 
are both devoid of hyperstaticity (force indeterminacy) 
and of hypostaticity (displacement indeterminacy), apart 
from possible trivial motions in which all force-carrying 
grains move as one rigid body. These two properties have 
different origins [32], and are not valid under the same 
assumptions. The absence of hyperstaticity (h = with 
the notations of Sec. Ill C[) results from the generic disor- 
der of the packing geometry. It would hold true for arbi- 
trarily shaped rigid particles interacting by purely nor- 
mal contact forces whatever the sign of those forces, and 
it applies to the whole packing, whatever the contacts 
the rattlers might accidentally have with the backbone. 
The absence of hypostaticity property (except for trivial 
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mechanisms, k = fco), on the other hand, is only guar- 
anteed for spherical particles with compressive forces in 
the contacts, and it applies to the sole backbone. 

Due to the isostaticity property, the coordination num- 
ber should be equal to 6 in the rigid limit on the back- 
bone. If Nc is the number of force- carrying contacts, 
then the global (mechanical) coordination number is 
2Nc 

z — (possible contacts of the rattlers are discarded) , 

n 

and the backbone coordination number is defined as 

2Nr z 

than z, has the 



coordination 

z 

z , rather 



n(l-xo) l-xo 
limit 6 as K ^ +oo. In A samples (/t = 39000) one 
has z* = 6.074 ± 0.002 (and hence z ~ 5.995), the excess 
over the limit z* — 6 resulting from contacts that should 
open on further decreasing the pressure. 

The isostaticity property can be used to evaluate the 
density increase due to finite particle stiffness. To first 
order in the small displacements between P = (or 
K = -\-oo ) and the current finite pressure state A, one 
might use the theorem of virtual work js^l , with the dis- 
placements that bring all overlaps /i^ to zero, and the 
current contact forces. Such motions leading to a simul- 
taneous opening (/ly = 0) of all contacts are only possi- 
ble on networks with no hyperstaticity, because there is 
no cornpatibility condition on relative normal displace- 
ments [33 . This yields an estimate of the increase of the 
solid fraction A<I> over its value $0 in the rigid limit, as 



1 v-^ A* 



(29) 



This equality can be rearranged using the Hertz contact 
law ID) to relate Nij to htj, and relation (22). We denote 
as Z{a) the moment of order a of the distribution of 
normal forces Nij, normalized by the average over all 
contacts: 



Z{a) = 



(iV") 
(TV)"' 



(P5)) can be rewritten as: 



A$ = 35/3$i/3^(5/3) r 



7r\2/3 



(30) 



(31) 



In the isostatic limit which is approached at large n, the 
force distribution and its moments are determined by 
the network geometry, and we observed Z(5/3) = 1.284. 
Taking for z and $ the values at the highest studied 
stiffness level k {k — 39000), this enables us to evalu- 
ate the density change between those configurations and 
the rigid limit as A<I> ~ 1.15 x 10^*. As announced be- 
fore this is smaller than the statistical uncertainty on 
$, and hence this does not improve our estimation of 
the solid fraction <I>o of the packing of rigid particles 
(k = -l-cxo). Recalling = (P/E)'^/^, ([31]) means that 
the macroscopic relation between pressure and density 
has the same power law form (P cx (A<I>)^/^) as the con- 
tact law {N (X h^/"^). This was observed by OSLN. It 



would hold, because of the isostaticity property in the 
rigid limit, for whatever exponent m in the contact law, 
the prefactor of the macroscopic relation P cx (AC')'" 
involving Z{\-\- 1/m), a moment of the geometrically de- 
termined force distribution. 



As a consequence of ([3T|) . one can simply relate the 
bulk modulus of frictionless packings to the pressure, as 
observed by OSLN too, a property which will be used and 
discussed in paper III [3a|, which deals with elasticity of 
packings. 



Force distribution 



The force distribution we observe in A samples at high 
K values approaches the one of a rigid packing, which 
due to isostaticity is a purely geometric property. It is 
represented on Figure [21 The data presented here are av- 
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FIG. 2: (Color online) Probability distribution function P{f) 
of normalized contact forces / = N/{N} in A and A' configu- 
rations at high K. The dashed line is the fit proposed by DTS: 
P(/) = (3.43/^ -\- 1.45 - 1.18/(1 + 4.71/)) exp(-2.25/). 



eraged over 5 samples. Because of relation ([22|) . all sam- 
ples prepared at the same pressure have the same aver- 
age force, and this restores the "self-averaging" property, 
which OSLN observed to be lacking on using solid frac- 
tion instead of pressure as the control parameter. The 
choice of $ as a state variable, because of the finite size 
of the sample causing fluctuations of the threshold $0 
where P vanishes, is less convenient in that respect. 

Fig. [5] also shows that the form proposed by DTS to 
fit their data is in very good agreement with our results, 
except perhaps for large forces, for which it is a better fit 
for A' data - thus providing additional evidence that A' 
samples are closer than A ones to the DTS results. 
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F. Geometric characterization 

1. Pair correlation function 

Pair correlations should preferably be measured either 
with method I or method II, as there is no reason to elimi- 
nate rattlers before studying geometric properties. Com- 
parisons between pair correlation functions (r) and 
g^^{r) (Fig. [3]) show very little difference on the scale of 
one particle diameter. Results of Fig. [3] are very similar 
to other published ones {e.g., in DTS), with an apparent 
divergence as r ^ a and a split second peak, with sharp 
maxima at r = aV3 and r — 2a. The pair correlation 




FIG. 3: Pair correlation functions (r) and g'^ (r) versus r /a 
in A samples at k = 39000. Both definitions coincide on this 



scale (only the peak for r 



is slightly different). 



function should contain a Dirac mass at r = a in the 
limit of rigid grains, which broadens into a sharp peak 
for finite contact stiffness. The weight of this Dirac term 
or peak in the neighbor intercenter distance probability 
distribution function 47r|^r^(7(r) — 2A{r^ / a?)^g{r) is co- 
ordination number 2;, and the shape of the left shoulder 
of the peak at finite k is directly related to the force 
distribution P{N): 



(For (5 > 0) g{a -5) = 



48$(a- (5)^ 



This explains the observation by OSLN [3l| of the width 
of the g{r) peak decreasing approximatively as A(E>, while 
its height increases as (A$)~^, as the threshold density 
$0 is approached from above. The form of the distribu- 
tion of contact forces, which is determined by the geom- 
etry of the isostatic backbone, remains exactly the same 
for all small enough values of A<I>, with a scale factor 
proportional to A$^/^ , due to Eqns. and El] 

The sharp drop of g{r) at r = a-\/3 and r = 2a was 
found by DTS to go to a discontinuity in the rigid particle 
limit. This can be understood as follows. Each sphere 
has a number of first contact neighbors [z on average) 



at distance r = a if the grains are rigid, and a number 
of second contact neighbors (i.e. particles not in contact 
with it, but having a contact with at least one of its 
first contact neighbors). Such second contact neighbors 
will make up for a significant fraction of particles with 
their centers at a distance r < 2a, but none of them can 
be farther away. Futhermore, this leads to a systematic 
depletion of the corona 2a < r < 2a + 6 (with < 6 < a) 
by steric exclusion. 



As 



2. Near neighbor correlations 



pair correlations are conveniently ex- 



pressed with the gap-dependent coordination number 
z{h). z{h) is the average number of neighbors of one 
sphere separated by an interstice narrower than h. z(0) is 
the usual contact coordination number z. Function z(h) 
has three possible different definitions z^ , z^' and z^^^ 
according to the treatment of rattlers. All three of them 
were observed to grow as z(0) -I- ChF''^ for h smaller than 
about 0.3a, constant C taking slightly different values for 
z^, z^^ and z^^^ . z^^^{h) is equal to z* for h = 0, and is 
very well fitted with the value C = 11 found by DTS (s^ . 
Fig. 8]. z(h) deviates from this power law dependence 
corresponding to the rigid limit for small h, of the or- 
der of the typical overlap shown on Fig. 31 This 
power law corresponds to g{r) diverging as (r — a)""'^ 
as r — > a"*". Silbert et al., in a recent numerical study 



10 



— z (h), K=39000 
--- z"(h),K=18000 

■ z"(h), K=8400 

— slope 0.6 




FIG. 4: (Color online) Coordination number of near neighbors 
function of interstice h. The power law regime extends to 
smaller and smaller h values as k increases. 

of states with high levels of rigidity [79;] {n > 10^), re- 
port observing z^{h) to grow with an exponent closer to 
0.5, although somewhat dependent on the choice of the 
h interval for the fit. However this does not contradict 
our main conclusion that different numerical approaches 
track the same RCP state in the rigid limit. 
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3. Other properties of contact networks 

Ignoring rattlers (method III) one may record the den- 
sity of specific particle arrangements in the backbone. 
We thus find the contact network (joining all centers of 
interacting particles by an edge) to comprise a number 
of equilateral triangles, such that on average each back- 
bone grain belongs to 2.04 ± 0.04 triangles. In the rigid 
limit this gives a Dirac term for = 7r/3 in the distri- 
bution of angles 9 between pairs of contacts of the same 
grain. Tetrahedra are however very scarce (as observed 
by DTS), involving about 2.5% of the beads, and pairs 
of tetrahedra with a common triangle are exceptional (5 
such pairs in 5 samples of 4000 beads). Pairs of triangles 
sharing a common base are present with a finite density, 
which explains the discontinuous drop at r = aVS of 
g{r), this being the largest possible distance for such a 
population of neighbor pairs. 



G. Conclusions 

We summarize here the essential results of Section Hill 
about frictionless packings. 



1. Uniqueness of the RCP state 

Our numerical evidence makes a strong case in favor 
of the uniqueness of the simulated rigid packing state 
made with frictionless spheres under isotropic pressure. 
Specifically, we observed quantitative agreements with 
other published results |31i . j)2j in the coordination num- 
bers, the force distributions, the pair correlations and 
the frequency of occurrence of local contact patterns, 
even though different numerical methods have been used. 
The small remaining differences in solid fraction, propor- 
tion of rattlers, and probability of large contact forces all 
correlate with the duration of agitated assembling stage, 
which can be measured in terms of numbers of collisions 
per grain at a given density. This duration directly cor- 
relates to the packing fraction and to the small amount 
of crystalline order in the samples. We therefore checked 
in an accurate, quantitative way the traditional views 
about random close packing (RCP). The RCP state can 
be defined in practice as the unique state in which rigid 
frictionless spherical beads assemble in a static equilib- 
rium state under isotropic pressure, in the limit of fast 
compression, so that the slow evolution towards crystal- 
lization has a negligible influence. The Lubachevsky- 
Stillinger algorithm tends to produce packings with a 
small but notable crystalline fraction. 



2. Relevance of MD simulations, role of micromechanical 
parameters 

The uniqueness of the RCP implies that dynamical pa- 
rameters C, and / have no influence on the frictionless 
packing structures, at least in the limit of fast compres- 
sion rates. 

Standard MD methods compare well with specifically 
designed methods that deal with rigid particles, and 
prove able to approach the rigid limit with satisfactory, if 
admittedly smaller, accuracy. Recalling that n — 39000 
corresponds to glass beads under 10 kPa, it seems that 
laboratory samples under usual conditions might in prin- 
ciple (if friction mobilization can be suppressed) ap- 
proach the ideal (rigid particle) RCP state. 

Moreover, the time scales to assemble samples in MD 
simulations, if compared to estimated preparation times 
in the laboratory with such techniques as controlled plu- 
viation, has the right order of magnitude. This means 
that the assembling proces is rather fast in experimental 
practice when grains are deposited under gravity, which 
explains why densities above RCP are not directly ob- 
tained. Of course, in practice, many procedures produce 
anisotropic states. Anisotropic packings of rigid, fric- 
tionless balls, under other confining stresses than a hy- 
drostatic pressure, should differ from the RCP state, and 
the numerical simulations of gravity deposited packings 
of frictionless beads of refs. [2J, [2^ could be analysed in 
this respect. We chose here to study ideal preparation 
methods, and we only deal with isotropic systems. 

3. Approach to isostaticity in the rigid limit 

We checked that bead packings under typical labora- 
tory pressures such as 10 kPa might closely approach the 
isostaticity property of rigid frictionless packings. We 
showed that some observations made by O'Hern et al. [3l| 
on pressure or bulk modulus dependence on density, and 
on the shape of the first peak of g{r) were direct conse- 
quences of this remarkable property. 

IV. LOW-PRESSURE STATES OF FRICTIONAL 
PACKINGS OBTAINED BY DIFFERENT 
PROCEDURES 

A. Introduction 

It is well known that the introduction of friction in 
granular packings tends in general to reduce density and 
coordination number, as observed in many recent numer- 
ical simulations (see e.g. [13, 25, 26, 27]), and that fric- 
tional granular assemblies, unlike frictionless ones, can 
be prepared in quite a large variety of different states. 
In the field of soil mechanics, sand samples are tradition- 
nally classified by their density [ll|, , which deter- 
mines behaviors that have been observed in simulations 
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of model systems as well [2ll |27[. (Inherent anisotropy 
of the fabric, i.,e., the one due to the assembling pro- 
cess, rather than induced by anisotropic stresses, is a 
secondary, less influential state variable [l3j[3])- Engi- 
neering studies on sands usually resort to a conventional 
definition of minimum and maximum densities, based on 
standardized procedures [80| . 

The motivation of the present study is to explore the 
range of accessible packing states, as obtained by dif- 
ferent numerical procedures that produce homogeneous 
and isotropic periodic samples. We therefore chose to by- 
pass the painstaking computations needed to mimic ac- 
tual laboratory assembling methods, but we argue that 
our procedures produce plausible structures with similar 
properties. 

One key result is that density alone does not determine 
the internal state of an isotropic packing, because the 
coordination number can vary independently. 

The A-type configurations obtained without friction in 
Section IIVBI are local density maxima in configuration 
space (see Sec. IIII ClT) . Hence compaction methods can 
be regarded as strategies to circumvent the mobilization 
of intergranular friction forces. Two such procedures are 
studied here, in a simplified, idealized form: lubrication 
and vibration. We also simulate, as a reference, a state 
which can be regarded as a loose packing limit, at least 
with a definition relative to one assembling method and 
friction coefficient /i = 0.3 ; and we prepare, as an in- 
teresting limit from a theoretical point of view, infinite 
friction samples. 

Assembling procedures are described in Section flVBi 
geometric aspects are studied in Section flV CI and con- 
tact network properties in Section IIVDI Section IIV El 
summarizes the results. 



B. Assembling processes for frictional grains 

Just like in the frictionless case, for each one of the 
packing states, we prepare 5 samples of 4000 beads, over 
which results are averaged, error bars corresponding to 
sample-to-sample fluctuations. The equilibrium criteria 
are those of Section [ill B 21 supplemented with a similar 
condition on moments. To identify rattlers, we use the 
procedure defined in Section fill E 11 which is adapted to 
the case of frictional grains: spheres with as few as two 
contacts may carry forces (even large ones, as we shall 
see) and should not be regarded as rattlers. 

1. Looser packings compressed with final friction coefficient 

We used direct compression of a granular gas in the 
presence of friction (u = 0.3), another standard numer- 
ical procedure [H, [IJ, HI], in which the obtained den- 
sity and coordination number are decreasing functions 
of jjL [H, m, H^, [131. This produces rather loose sam- 
ples hereafter referred to as D (B and C ones, to be 



defined further, denoting denser ones, closer to A, but 
arguably more "realistic"). D samples were made with 
exactly the same method as A ones (see Sec. IIII B 2p . 
except that the friction cocfhcicnt ~ 0.3 was used in- 
stead of /i = 0. In principle, D configurations should 
depend on initial compaction dynamics: increasing the 
rate of compression could produce denser equilibrated 
packings, just like a larger height of free fall, whence 
a larger initial kinetic energy, increases the density of 
configurations obtained by rain deposition under grav- 
ity [67]. We request the reduced compression rate /, de- 
fined in (|23p , not to exceed a prescribed maximum value 
-fmax- The choice of /max = 10~^ and C, = 0.98 yields 
solid fractions $ = 0.5923 ±0.0006 and backbone coordi- 
nation numbers z* = 4.546±0.009, with a rattler fraction 
xo = (11.1±0.4)%. These data correspond to P = 1 kPa 
(or K ~ 181000). Very similar results are obtained on 
using a different, but low enough pressure, such as 10 
or even 100 kPa, as remarked in [4l| (where 2D sam- 
ples were assembled by oedometric compression), and as 
indicated in table IIIIl However, a quasistatic compres- 
sion from P ~ 1 kPa to 10 or 100 kPa produces slightly 
different states at the same pressure. The influence of C 
should disappear in the limit of slow compression, / — s- 0. 
A practical definition of a (/.t-dependent) limit of loose 
packing obtained by direct compression, can therefore 
be proposed as the / ^ limit of our D states. As 
reported in table IIIIl a value of the damping parameter 
ten times as small as the standard one C = 0.98 results in 
quite similar configuration properties, on compressing a 
loose granular gas under P = 10 kPa, with /max — 10^^. 
So did in fact faster compressions, with /max = 10"^ , 
keeping C — 0.98. The data of Table IIIIl thus suggest 
that we very nearly achieved the independence on dy- 
namical parameters that is expected in the I —^ limit 
with our choice of control parameters. We note, however, 
that other possible definitions of a random loose packing, 
such as the one by Onoda and Liniger [81i] result in dif- 
ferent (smaller) solid fractions. Looser arrangements of 
equal-sized spherical particles can also be stabilized with 
adhesive contact forces, e.g. on introducing the capillary 
attractions produced by the menisci formed by a wetting 
fluid in the interstices between neighboring grains [82!|. 

In addition to packing fraction $, coordination num- 
ber z*, fraction of rattlers xq, Tabic [Till lists the reduced 
second moment Z{2) of the normal force distribution, as 
defined in pO)) , the proportion of two-coordinated beads 
(to be discussed in Section lTVDI) . X2, and the average val- 
ues of ratios ||T||/A^ (friction mobilization) among con- 
tacts carrying normal forces larger and smaller than the 
average, respectively denoted as Mi and A/2. As a result 
of some amount of quasistatic compression of the initial 
assembly, the width of the force distribution decreases, 
as witnessed by smaller values of Z{2) in table [Till and 
so does the mobilization of friction, as measured by Mi 
and A/2. The effects of compression on the structure and 
the forces are further studied in paper II [37| . On com- 
paring numerically simulated loose packings to experi- 
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TABLE III: Isotropic states of type D, from direct compression of the granular gas at the indicated pressure (rows marked 
"gas"), or from gradual, quasistatic compression (rows marked "QS") of solid samples made at the lowest pressure 1 kPa (or 
K ~ 181000). Tests of the influence of viscous dissipation parameter ( and maximum value /max of reduced strain rate in 
compression are also made for configurations compressed from a granular gas to 10 kPa. 



Origin 


P (kPa) 


c 


^max 




z* 


xo (%) 


X2 (%) Z{2) 


Ml 


A/2 


gas 


1 


0.98 


10" 


3 


0.5930 ± 0.0007 


4.546 ± 0.009 


11.1 ±0.4 


2.39 


1.58 


0.160 


0.217 


gas 


10 


0.98 


10" 


3 


0.5946 ± 0.0006 


4.59 ± 0.02 


10.2 ± 0.2 


2.07 


1.59 


0.159 


0.213 


gas 


10 


0.098 


10" 


3 


0.5938 ± 0.0008 


4.61 ± 0.02 


10.9 ±0.2 


1.79 


1.57 


0.150 


0.194 


gas 


10 


0.98 


10" 


-1 


0.5931 ± 0.0002 


4.60 ±0.01 


10.2 ± 0.7 


1.80 


1.59 


0.159 


0.212 


QS 


10 


0.98 


10" 


3 


0.5931 ± 0.0006 


4.641 ±0.011 


10.1 ±0.4 


2.33 


1.46 


0.146 


0.188 


gas 


100 


0.98 


10" 


-3 


0.5975 ± 0.001 


4.69 ± 0.02 


8.9 ±0.5 


1.66 


1.61 


0.153 


0.197 


QS 


100 


0.98 


10" 


-3 


0.5936 ± 0.0006 


4.79 ± 0.02 


8.6 ±0.4 


2.05 


1.40 


0.138 


0.178 



mcnts, it should be recalled that samples are assembled 
under low pressures in the laboratory: the hydrostatic 
pressure under a 1 cm thick layer of glass beads is about 
0.15 kPa. Numerical configurations under higher confin- 
ing pressures corresponding to mechanical tests in the 
laboratory {e.g., sound propagation) are more appropri- 
ate models if the testing pressure is significantly larger 
than the initial, assembling pressure - as for the "QS" 
samples of table IIIII 

The effects of such proportions of rattlers in granular 
packings as reported in table IIIII have to our knowledge 
never been studied in detail. It should be emphasized 
that this relatively large population of rattlers does not 
jeopardize the global stability of equilibrium configura- 
tions, as the stiffness matrix of the force-carrying network 
is found devoid of floppy modes (apart from harmless, lo- 
calized ones associated with two-coordinated particles, to 
be discussed in Section IIVD[) . 

Our D samples should be compared with the simula- 
tions reported by Zhang and Makse [SS^, in which loose 
sphere packings were also prepared by isotropic compres- 
sion. Those authors observed, in some cases, lower pack- 
ing fractions than D values, $ ~ 0.57. Their assembling 
method is however different: they use a strain-controlled 
procedure, with a constant compression rate, and then re- 
lax the final state at constant volume. In this approach, 
the pressure reaches very high levels, several orders of 
magnitude as large as the final value, before samples fi- 
nally stabilize [83, Fig. 3]. Zhang and Makse report 
some dependence of the final state on the compressing 
rate. Once translated into the dimensionless parameter 
I we have been using here, strain rates used in |83| . de- 
fined with the typical pressure value P = 100 kPa, range 
between / = 0.1 and I = 100. The slowest compres- 
sion reported in [s^ is therefore 100 times as fast as the 
upper limit for e we have been enforcing in this work. 
Viscous forces also differ between the present simulations 
and those of Ref. [s^ , in which "global damping" terms 
are used (i.e., forces opposing the individual motion of 
particles, rather than relative motions). 



2. Use of low friction coefficients: imperfect lubrication 

One way to limit the effects of friction consist in lubri- 
cating the grains, as in the experimental study reported 
in [8J|. If all intergranular friction could be suppressed 
in the assembling procedure, i.e., for perfect lubrication, 
then the structure of isotropic packings would be the one 
denoted as A, studied in Section ilTTl The effect of a small 
friction coefficient in the contacts while the grain assem- 
bly is being compressed can be regarded as a crude, sim- 
plified model for imperfect lubrication. We made sam- 
ples, denoted as B, by compressing the granular gas, just 
like in the A and D cases, with /x = 0.02. In order to 
approach the limit of slow compression rates better, we 
started from D configurations, decreased the friction co- 
efficient to ^ = 0.02, and then requested that / < 10"''^ 
while the samples got further compressed to equilibrium 
under 1 kPa (k ~ 181000). (In view of the results in the D 
case of Section llVB II we do not expect the final B state 
to be sensitive to damping parameter (.) We observed 
that this small friction coefficient had a notable effect on 
the final solid fraction, as the value $ = 0.6270 ± 2.10""* 
is significantly below the frictionless (A) result, while the 
coordination number on the active structure is slightly 
reduced, down to z* — 5.75 ± 8.10""^, and the fraction of 
rattlers raised sHghtly, to xq = (1.95 ± 0.02)%. 



3. Dense, fractional packings obtained by shaking 

Another practical strategy to obtain dense configura- 
tions is to shake, vibrate or apply repeated "taps" on 
granular samples Issl [3^ . Such procedures involve the 
introduction of kinetic energy into already quite dense 
assemblies. In order to investigate their possible effects 
at a limited computational cost, we avoid the direct sim- 
ulation of repeated shakes and adopted the following pro- 
cedure. Starting from the dense A configurations (made 
without friction and described in Section IIII[) , we first 
apply a homogeneous expansion, multiplying all coordi- 
nates by a common factor A slightly larger than 1 . With 
equilibrated A states under confinement level k — 39000, 
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the chosen value A = 1.005 is more than enough to sep- 
arate all pairs in contact. Then, in order to mimic, in 
an idealized way, the motion set up by a shaking ex- 
citation, the beads are given random velocities (chosen 
according to a Maxwell distribution) , and interact in col- 
lisions which preserve kinetic energy, while the volume of 
the cell is kept constant. This "mixing" stage is simu- 
lated with the "hard sphere molecular dynamics" (event- 
driven) scheme (just like our initial granular fluids are 
prepared at $ = 0.45, as described in Section UlIB 2p . 
It is pursued until each particle has had ricou = 50 col- 
lisions on average. The final preparation stage is a fast 
compression: velocities are set to zero, particles regain 
their elastic and dissipative properties (as defined in Sec- 
tion III A| with friction coefficient /i = 0.3, and viscous 
dissipation, ^ = 0.98), the external pressure P = 10 kPa 
is applied via the deformable periodic cell, until a final 
equilibrium is reached. 

The final state is hereafter referred to as C. Quite un- 
surprisingly, its solid fraction, $ = 0.635 ± 0.002, stays 
very close to the RCP value obtained in the A state. 
However, the coordination number is considerably lower, 
z* = 4.56 ± 0.03, which is as small as the value obtained 
in the loose ($ ~ 0.593) D state, while the proportion 
of rattlers raises to xq = (13.3 ± 0.5)%. Remarkably, on 
comparing states B and C, the latter has a higher den- 
sity, but a much lower coordination number, and a much 
higher fraction of rattlers. 

We did not thoroughly investigate the infiuence of pa- 
rameters A and ncoii, introduced in the preparation proce- 
dure, on the resulting C states. In the following we focus 
on configurations obtained with the values A = 1.005 and 
ricoU = 50. Yet, we noted that an increase of A to 1.01 
entailed only very slight changes of and z* (which re- 
spectively decreased to 0.633 and 4.54), and that the sup- 
pression of the "mixing" stage (i.e., setting ncou to zero) 
resulted in much higher z* values (around 5.5). Like- 
wise, we did not check for a possible effect of C on the 
final state. Smaller values than the large one C = 0.98 
used in our simulations of the final compression stage of C 
sample preparation are likely to have analogous effects to 
an increase of the duration of the agitated mixing stage, 
and should not increase the final coordination number. 



4- Global state variables: summary and discussion 

Table IIVI gathers some of the parameters characteriz- 
ing the four different packing states studied in the present 
paper. (Mi and M2 were defined in connection with ta- 
ble |TTT]). From table HVl configurations with low coordi- 
nation numbers (C and D) appear to exhibit a somewhat 
wider normal force distribution (as measured by Z{2)), 
and a significant mobilization of friction, with typical val- 
ues of ||T||/7V around = 0.15 for larger than aver- 
age normal force components N (Mi), and significantly 
above ^/2 for smaller N values (M2). 

The existence of C states shows that there is no system- 



atic relationship between density and coordination num- 
ber, contrarily to some statements in the literature [7|. Of 
course, for one particular assembling method both quan- 
tities will often vary in the same direction as functions of 
some control parameter. For instance, on preparing sam- 
ples by deposition under gravity, both density and coor- 
dination number are increasing functions of the height 
of free fall [67]. However, our results show that different 
preparation methods might lead to contrasting results. 

Our results about density and coordination number 
can be likened to observations made before in numeri- 
cal models of sphere packings obtained with geometric 
construction rules [s^ . The simplest versions of such al- 
gorithms m, , which mimic deposition under gravity, 
add particles one by one by dropping and rolling them 
in contact with one or two previously deposited parti- 
cles, until they are fixed when they rely on three con- 
tacts. Those produce packings with z = 6. More re- 
fined versions thereof 85, 88J also involve other, more 
collective types of moves. The final configurations then 
contain "bridges" or "arches" [s^ , defined as sets of par- 
ticles the final stabilization of which is mutual and collec- 
tive. In such arches, each grain relies on three others, but 
some pairs mutually rely, in part, on each other. Those 
"bridged structures" have much lower coordination num- 
bers, down to about 4.5. 

It is not clear, though, to what extent our config- 
urations, which were obtained within a full mechani- 
cal model, compare to those that result from such ap- 
proaches. As shown, e.g., in [90], deposition algorithms 
based on geometrical rules are supposed to ensure lo- 
cal stability properties, but the resulting granular pilings 
might turn out to be globally unstable. Moreover, a de- 
scription of our packings as a sequence of arches placed 
one after another, assuming it is conceivable, would seem 
to contradict their homogeneity and isotropy: it is rather 
arbitrary, in isotropic packings, to regard some parti- 
cles as "relying" on some others. Such a description was 
therefore not attempted. 

C. Geometric characterization 

1. Pair correlation functions 

As observed in previous experimental [tI] and numer- 
ical [l^] results, pair correlation functions present the 
same features at lower densities as at the largest one 
<I> ~ 0.64, in a weakened form, as shown on Fig. [51 On 
comparing those functions for states A to D, we observed 
what follows. 

• C samples, obtained from A ones after small rear- 
rangements, exhibit pair correlations that only dif- 
fer in the detailed shape of the peaks {e.g. below 
1.05a), and is indistinguishable elsewhere. 

• In spite of the large number of rattlers in samples C 
and D, {r) and g^^{r) (as defined in Section llll E|) 
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TABLE IV: Isotropic states (k ~ 39000 for A and C, k ~ 181000 for B and D) for different assembling procedures. 



Procedure 


$ z* xo (%) X2 (%) Z(2) Ml Ma 


A 


0.6370 ± 0.0002 6.074 ± 0.0015 1.3 ±0.2 1.53 


B (fio = 0.02) 


0.6271 ± 0.0002 5.80 ± 0.007 1.95 ± 0.02 ~ 10"'' 1.52 0.016 0.018 


C (A = 1.005) 


0.635 ±0.002 4.56 ±0.03 13.3 ± 0.5 2.64 1.65 0.135 0.181 


D 


0.5923 ± 0.0006 4.546 ±0.009 11.1 ± 0.4 2.39 1.58 0.160 0.217 



1 1.2 1.4 1.6 1 




FIG. 5: Pair correlation functions g{r) (definitions (r) and 
g^'{r) coincide on this scale) for configurations B, C, D. 



cannot be distinguished on the scale of Fig. \5\ 



• The depth of the trough between r/a = 1.1 and 
r/a = 1.5 increases with density. 



2. Near neighbor coordination numbers 

The gap-dependent coordination number z^^{h) is 
shown on Fig. [6] for samples A to D. Fig. [6] shows that, 
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FIG. 6: (Color online) Coordination number for neighbors at 
distance < /i, z^'' (h), for configurations A (red, upper dashed 
line), B (blue, middle dashed line), C (black, solid line) and 
D (green, bottom dashed line). 

as might have been intuitively expected, z{h) correlates 
with coordination number for small h and with density 
for larger distances, h > 0.04a. We preferably use defi- 
nition z^^ , which is obtained on bringing the rattlers in 
contact with the backbone with small, random forces, as 
explained in Sec. lIIIEl z^^{h) can be thought of as more 
physically meaningful than z^{h), which directly results 
from the simulation of the packing, and is somewhat am- 
biguously defined because the positions of the rattlers are 
not specified. Functions z^^{h) corresponding to B and 
C states cross each other for h ~ 0.02a. 

z^ {h) and z^^{h) might be fitted by power laws: 



z'{h) = B' + Aih^' 



(32) 



• The integral below the peaks correlates with den- 
sity, but the height and sharpness of the drop at 
r/a = a/3 and r/a = 2 correlate with coordina- 
tion number (which is larger for B than for C), in 
agreement with the interpretation of such features 
suggested in Sec. IIIIFl 



Figs m [5] and El display z^{h) - B' and z"{h) ~ B 
as functions of h on logarithmic plots for samples D, B 
and C (due to the influence, at short distance, of contact 
deflections on z{h) data, flt parameters B^ and B^^ are 
a little smaller than z and ^■'^^(0)). For h values smaller 
than 10~^, the lowest limit on the axis on Figs. [71 El and [51 
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the gap h is of the same order as elastic deflections, k~^, 
and we do not observe simple power laws (as on Fig. [?]). 
These figures (on which the corresponding values of 
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0.1 - 



0.0001 




FIG. 7: (Color online) z\h) - z (black) and z^\K) ~ z"{0) 
(blue) versus h on doubly logarithmic plot for D samples 
at lowest pressure. The slopes of the corresponding dotted 
straight lines (power law fits) are — 0.51 and Pi = 0.35 
(see Eqn.|32}. 




0.0001 



FIG. 8: (Color online) Same as Fig. [71 for B samples. Dotted 
lines have slopes (3ii = 0.56 and /3/ = 0.45. 



.JI 



z^{0) and z^^(O) are also provided) show that and 
have quite different h dependences. This should be 
accounted for on studying the closing of contacts due to 
compression (see paper II 3J] ) . As a result of the com- 
putation leading to the equilibrated state, many pairs of 
neighbors end up separated by a very small interstice, 
so that z^ — z already reaches values larger than 0.3 for 
h — lO^^a in samples C and D. z^ then grows with h 
more slowly than z^^, with /3// > /?/, and a power law fit 
of lesser quality. Exponent /J//, which should be regarded 
as a more intrinsic quantity than /?/, appears to correlate 
with solid fraction. It has the same value 0.6 in samples 




0.0001 



FIG. 9: (Color online) Same as Fig. [T] for C samples. Slopes 
of dotted straight lines: /3// = 0.6 and (3i — 0.37. 



C and A (Figs [9] and |4]), decreases to 0.55 in the inter- 
mediate density state B, and to 0.51 for the least dense 
one, D. The power-law form of z^^ extends to h ~ 0.3 in 
configuration A, to about h ~ 0.2 in configuration B, and 
only to 0.04 and 0.05 for C and D. Such limited power 
law ranges preclude comparisons with the experimental 
data of |3l. 

One has z^(0) ~ z as a very good approximation, since 
contacts carrying no force in the equilibrated state ob- 
tained by MD are very scarce. z^^{0), on the other hand, 
is the geometric coordination number once all rattlers are 
pushed against the backbone. It is larger than the me- 
chanical coordination number, z. Specifically, because 
all rattlers, in treatment II, are dealt with as frictionless, 
one has: 



(0) = 



6xn - 



2Nr. 



(33) 



in which Nrr is the final number of contacts between 
rattlers once they are positioned against the backbone. 
Value of z-'^^(O) can be read on Fig. [HI As pointed out 
m Section ImETI they can be compared to coordination 
number values of samples under gravity. The results of 
Silbert et al. [2J, Figs. 2 and 3], with jj, = 0.3, corre- 
spond approximately with our D samples: $ ~ 0.59 and 
2:^^(0) ~ 5. The systems simulated in [IJ] are however 
not isotropic. None of the samples made under gravity 
in [23, appears to resemble our C state. 



3. Absence of crystalline order, local order parameters 

All indicators of incipient crystalline order given in Ta- 
ble |T] for frictionless A samples take lower values in states 
B and D, while C configurations, due to their geometric 
proximity, are close to A ones in this respect. Already 
scant in dense configurations assembled without friction, 
traces of crystallization are thus negligible in looser ones 
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obtained with frictional beads. Like for A samples, nu- 
merical data on configurations around one sphere i, as 
characterized by the pair (Q4(«), QeC*)); used by Aste et 
al. Q, are presented for states B, C, and D and com- 
pared to their experimental results in Appendix |DJ It 
is observed that local disordered environments around 
one grain are very similar in numerical and experimen- 
tal configurations of equal densities, while local HCP-like 
arrangements occur with similar (low) frequencies, and 
FCC-like ones are present in the laboratory, but not de- 
tected in simulations. 



D. Properties of force networks 

1. Local contact coordination numbers 



The distribution of local coordinations is given in Ta- 
ble |V1 for both mechanical (I) and geometric (II) defini- 



State 


Xo Xl X2 X-i Xa X5 X(i X7 Xs Xg XU) 


B (I) 


1.95 0.05 0.5 16.3 26.9 27.5 18.4 7.1 1.4 0.1 


C (I) 


13.3 2.6 15.1 26.5 23.4 13.2 4.8 0.9 0.15 


D(I) 


11.1 2.4 13.8 29.1 25.6 13.3 4.0 0.7 0.03 


B (II) 


2.1 15.3 26.5 27.5 18.9 7.8 1.6 0.2 


C (II) 


18.9 22.8 27.0 18.6 9.3 2.8 0.5 0.1 


D(II) 


17.2 25.4 29.0 18.5 7.7 1.9 0.3 



TABLE V: Percentage Xi of grains having i contacts in con- 
figurations B, C and D, on ignoring contacts with or between 
rattlers (I), or on fixing them onto the backbone with small, 
random forces (II). 



tions of contacts, for states B to D. Compared to the A 
case (Table HI)), the distribution is shifted to lower values, 
with 4 and 5 the most frequent ones (rather than 6) in low 
coordinated packings C and D. Those samples also have 
quite a large population of three-coordinated spheres, 
and a notable one of divalent (two-coordinated) parti- 
cles. This contrasts with the frictionless case for which 
X2 = Xs = 0. Without friction, divalent spheres, from 
Eqn. [201 written with 7V^ = 2, iV/=3,/i = 0, would im- 
ply a mechanism and therefore an instability, and three- 
coordinated ones, in the absence of external forces and 
cohesion, cannot be equilibrated by non-vanishing nor- 
mal forces the net effect of which necessarily pushes them 
away from the plane defined by the three centers of their 
touching neighbors (the non-generic case with the four 
sphere centers within the same plane leading to an insta- 
bility). Spheres with three contacts therefore need some 
mobilization of friction to transmit non- vanishing forces 
in an equilibrium configuration, for tangential compo- 
nents are requested to cancel this net repulsion. The 
Coulomb condition then restricts such possible configu- 
rations to flat enough tetrahedra for contact forces to 
remain within the friction cone. This explains the small 
value of ^3 in low friction {fi = 0.02) B samples. 



2. The special case of divalent grains 

With friction, the small structure formed by one sphere 
having two contacts with fixed objects fFig. [TU)) . due to 
Eqn. 1191 in which the number of degres of freedom (6) 
is equal to the number of contact force coordinates, has 
a degree of force indeterminacy equal to its number of 
independent mechanisms: H = fc. In fact, both numbers 
are equal to I . Self-balanced contact forces (see first part 
of Fig. [ro]) are oriented along the line joining the two con- 
tacts, just like in the corresponding 2D case dealt with 
in [431, and their amplitude is a free parameter (the de- 
gree of "wedging" of the grain in the corner formed by 
its two neighbors [i^l). Such a possibility requires in 
practice that the angle, which we denote as a (Fig. [TU)) . 
between the line joining the centers of 1 and 2 and the 
one joining the contact points be smaller than the angle 
of friction, for the total contact force to stay within the 
Coulomb cone. In C and D samples, we observed tana to 
be distributed rather evenly between and /i, while the 
intensity of forces transmitted by divalent spheres ranged 
from to a few times Pa^. The mechanism associated 




FIG. 10: Equilibrium and free motion (mechanism) of one 
sphere (marked 1) with two contacts (with particles marked 
2 and 3). (a): in the plane of the three centers, normal and 
tangential components of the two contact forces balance along 
the line joining the contact points (dotted line), (b): seen 
from above, along the direction of the line of the centers of 
spheres 2 and 3, sphere 1 can move and occupy the diflferent 
positions depicted with dashed lines, its center describing the 
dotted circle around the 2-3 axis. 

with divalent spheres is a free rolling motion on the two 
contacts, the line joining the contact points being the 
instantaneous axis of rotation, as shown on fig. llOf b). 
In this motion, it is readily checked (see Appendix [E| 
that the rules given in Appendix [B] for the evolution of 
contact forces in the case of rolling and pivoting specify 
that contact forces will remain carried by the line join- 
ing the two contacts, with a constant intensity, as such 
contacts move on the surface of the fixed spheres. Such 
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forces will do no work and the kinetic energy of the mo- 
bile sphere as well as the elastic energy stored in its two 
contacts will be kept constant. The equilibrium of the 
divalent particle is thus marginally stable, with matrix 
K''^-' causing zero acceleration to the free motion. How- 
ever, such a motion does affect the balance of moments 
on spheres marked 2 and 3 on Fig. [TOl since the con- 
stant force is applied at a point that is moving on their 
surface. Therefore the stability of such free motions, as 
regards the global contact network, requires some addi- 
tional analysis - which is tackled in Appendix |E1 where 
it is concluded that the packing remains stable. Unlike 
in frictionless packings [s^, [MI ) mechanisms in the pres- 
ence of friction do not necessarily lead to instabilities. On 
building the constitutive stiffness matrix K^^^ in the sam- 
ples we studied, we could check that no other mechanism 
was present on the backbone than those rolling motions 
of divalent spheres. Once some stiffness element is in- 
troduced to impede the free motion of divalent spheres, 
one can e.g. check that the Cholesky factorization of the 
stiffness matrix only involves strictly positive terms on 
the diagonal. 

There can be no contact between two divalent spheres, 
as the simultaneous equilibrium of each of them with two 
forces carried by the line joining its two contact points, 
as on Fig. [TOl is impossible. 



3. Degree of force indeterminacy 

On the backbone, with n{l — a;o) spheres and = 
6n(l — xq) + 3 degrees of freedom, one has on average 
z*n{l — xo)/2 contacts and fc = 3 -I- X2n independent 
mechanisms, hence a degree of force indeterminacy (hy- 
perstaticity) , from p9p . given by: 



H = A^: 



with 



2X2 



3(1 -xo) 

2X2 



and 



(34) 



z + 



Hl-xo) 



The backbone is devoid of force indeterminacy (h = 0) 
when its coordination number is equal to Zq. Because 
of the mechanisms associated with divalent beads Zq is 
strictly smaller than 4. Alternatively, one can define a 
"corrected" backbone coordination number z**, as writ- 
ten in (f34l) . which is equal to 4 in the absence of force 
indeterminacy. As to the global mechanical coordination 
number zq corresponding to the absence of force indeter- 
minacy, its value, given by 



2X2 

zq = 4(1 - xo) —, 



(35) 



is well below 4. From the data of table|Vl zg is about 3.45 
in state C and 3.54 for D (while z is close to 4). Although 
H is relatively small compared to the number of degrees 



of freedom Nf, the samples with friction and low coordi- 
nation are still notably hyperstatic - a conclusion shared 
by other studies , which we reach here in the slightly 
different context of packings in a uniform state of stress. 
Unlike for frictionless sphere assemblies, there is actually 
no special reason to expect packings with intergranular 
friction to become isostatic in the rigid contact limit. The 
essential difference is that contact forces can no longer be 
regarded as enforcing hard geometric constraints like im- 
penetrability, and hyperstatic configurations do not re- 
quire exceptional arrangements or matching of particle 
sizes as in the frictionless case 32, tI, |76|, |92[. 

Unlike us, Zhang and Makse [s^ speculate that iso- 
static packings could be obtained in the limit of slow 
compressions of samples with ordinary values of /i. This 
is however due to a divergence of interpretation, rather 
than a contradiction in numerical results, since their min- 
imum coordination numbers z*, excluding rattlers, are 
similar to ours, z* ~ 4.5. Zhang and Makse could only 
approach configurations devoid of hyperstaticity on set- 
ting the friction coefficient to infinity (see Section HVD 61 
below). The degree of force indeterminacy per degree 
of freedom on the backbone is still equal, from ([M)) . to 
0.141 in D samples and 0.145 in C ones, and varies very 
little with compression rate in the range we explored, 
which extends to significantly smaller values than the 
ones used in [s^, as stressed above. It is not obvious 
whether special experimental situations might occur in 
which real granular assemblies approach vanishing de- 
grees of hyperstaticity. 



4. Distribution of normal forces 

Since the force-carrying structure maintains a non- 
vanishing degree of force indeterminacy even in the rigid 
limit for frictional packings, the force distribution in 
states B, C, and D, unlike in A configurations, is no 
longer a geometrically determined quantity in the rigid 
limit. 

The distribution of normal components of contact 
forces (normalized by its average (A^)) is shown on Fig. 1111 
for all 4 configurations A, B, C, and D, at the lowest 
pressure (as given in Table IIV[) , at the end of the as- 
sembling process. We observe, as in many other numeri- 
cal [25I, [73, m, HI] and experimental [9I, studies, an 
approximately exponential decay of P{f) for large values, 
which is somewhat slower in states with low coordination 
number (in agreement with the values of Z(2) given in 
Table HV)) . It should be pointed out, however, that much 
larger differences between force distributions in the four 
studied states A, B, C, D will appear on increasing the 
confining pressure (see [37[, paper II). This is already ap- 
parent in the dependence of Z{2) on the previous history 
of D samples in table IIIII 

All probability distribution functions show a local min- 
imum for / — > 0, except in state C. Although it was re- 
marked in past publications [97 | that an upturn of the 
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FIG. 11; (Color online) Distribution of normal forces, nor- 
malized by its average, / = -^j, in states A (red crosses), 
B (blue asterisks), C (black square dots) and D (green open 
squares) at the lowest simulated pressure. 



6. The limit of large frietion coefficients 

Motivated by the search for an isostatic hmit in pack- 
ings with intergranular friction, Zhang and Makse [s^ 
assembled numerical packings with friction coefficients 
equal to infinity. Then they could get z* ~ 4.15. In 
order to investigate, in paper III [sl] the elastic proper- 
ties of tenuous contact networks, we also prepared a set 
of 5 configurations on compressing a granular gas under 
P = 1 kPa with fi = -l-oo and condition / < 0.001. These 
states, herafter referred to as Z configurations, have solid 
fraction $ = 0.5917 ± 0.0008, backbone coordination 
number z* = 4.068 ± 0.006, proportion of rattlers xq — 
(18.4 ± 0.5)% and of divalent beads X2 = (6.8 ± 0.5)%, 

and H is indeed small in that case, -3^ — 0.031. It seems 

"f 

therefore very plausible that the degree of hyperstaticity 
vanishes in this case, as k — > oo, for slowly assembled 
packs. In the light of this observation, the absence of 
such a limit for finite /i can be attributed to sliding fric- 
tion destabilizing barely rigid structures, which collapse 
and tend to form contacts in excess over the minimum 
count. 



p.d.f. at low forces, as for C configurations, appeared 
when packings were not fully equilibrated, C configura- 
tions satisfy equilibrium requirements as well as the oth- 



E. Conclusions 



We summarize here the most salient results of Sec- 
tion |TVl about systems with friction. 



5. Friction mobilization 

As frictionless sphere packings are unstable for z* < 
6 (3^ . a certain level of friction mobilization is required 
in states B to D, even under isotropic stresses. In partic- 
ular, non- vanishing tangential forces are indispensable to 
ensure the equilibrium of grains with 2 and 3 contacts, 
for which they relate to the local geometric configuration, 
as discussed above (see e.g., Fig. [TO]). Some information 
about friction mobilization is given in Table IVII In Ta- 
ble IVII we distinguish between contacts carrying normal 
forces larger and smaller than the average. Those two 
populations of contacts, as first distinguished in [OTtlosj. 
as the "strong" and "weak" networks, are attributed dif- 
ferent roles, especially under anisotropic stresses. We 
merely use these categories here to gather information, 
in a compact, summarized form, about some aspects of 
force networks, which correlate to the force level. Ta- 
ble |VT] shows that, although three-coordinated particles 
tend to carry small forces on average, a notable propor- 
tion of them, and of the divalent ones, participate in the 
strong "force chains" with larger than average force lev- 
els. Friction mobilization is necessary in the contacts 
with such spheres, and reaches similar levels on average 
in the whole network. It it larger for contacts carrying 
small loads. 



1. Diversity of states of equilibrated packings 

The variety of inner structures of isotropic bead pack- 
ings we have obtained, as summed up in Table HVl shows 
that the solid fraction is not the only variable determin- 
ing the internal state of an isotropic equilibrated pack- 
ing. In particular, the backbone coordination number z* 
can vary throughout the whole interval from about 4.5 
to 6 in the rigid limit (k — > +00), in packings with a 
solid fraction close to the RCP value. Some systems can 
be denser than others, but with a considerably smaller 
coordination number. Systems compacted by vibration 
should have smaller coordination numbers than systems 
assembled with low friction coefficients. 



2. Geometry and length scales 

Geometric characteristics corresponding to length 
scales above about 4 or 5% of the particle diameter corre- 
late with density: this applies to the global shape and the 
area under the peaks of pair correlation function g{r), to 
near neighbor gap-dependent coordination number z{h), 
and to local neighbor arrangements around one bead, as 
measured by the local order parameters charted in Ap- 
pendix |D] 
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State 


Xi 


Ml 


M2 










Ar(3) xl^^ 


Mf^ 




B 


0.42 


0.016 


0.018 


n. d. 


n. d. 


n. d. 


n. d. 


0.23 0.09 


0.02 


0.01 


C 


0.41 


0.14 


0.18 


0.58 


0.19 


0.14 


0.17 


0.24 0.09 


0.13 


0.04 


D 


0.41 


0.16 


0.22 


0.71 


0.22 


0.16 


0.18 


0.26 0.10 


0.15 


0.05 



TABLE VI: For states B, C, and D, at the lowest pressure (see table ITV)l . proportion Xi of contacts carrying normal forces 
larger than the average (A''), average ratio for contacts with N > (N) (respectively, N < (N)) Mi (resp. M2). The same 
quantities with superscripts (2) and (3) apply to contacts implying spheres of coordination 2 and 3. N^^\ A''^' are the average 
normal forces carried by these two populations of contacts, normalized by Pa^. Note that divalent grains are absent in state 
B, and the corresponding quantities are therefore not defined. 



Features associated with smaller scales, such as z{h) 
for h < 0.04a or the shapes of the peaks of g{r) - part 
of which approach a discontinuity in the rigid limit - 
correlate with coordination number. 

Finally, a third, smaller length scale K~^a is associated 
with contact deflection and vanishes in the rigid limit. On 
this scale the geometric properties of the packing depend 
on the contact law. 



3. Role of rattlers 

Rattlers represent a significant fraction (above 10%) of 
the particles in poorly coordinated systems, although the 
packings are stable. The treatment of rattlers - whether 
they are left floating in arbitrary positions (I) or gen- 
tly pushed against the backbone by small forces (II) ~ 
changes geometric data on the intermediate scale men- 
tioned in the previous paragraph (Sec. IIVE 2p . such as 
the exponent of a power-law flt of z{h). Treatment II 
should be preferred if comparisons are to be made with 
packings under gravity. It leads to the definition of a 
"geometric" coordination number, z^^(O), different from 
the mechanical one (z). 

4-. Influence of micromechanical parameters 

Our data suggest (table IIIip that the states obtained 
on isotropically compressing a granular gas no longer de- 
pend on viscous damping parameter C in the limit of slow 
compression (/ — s- 0). If the true value of the friction 
coefficient is used at this stage (i.e., without "lubrica- 
tion"), the resulting state (our D configurations) can be 
regarded as a reference, loose packing limit of this as- 
sembling method. Other methods nevertheless result in 
lower densities. 



5. Force indeterminacy 

In samples assembled by isotropic compression, the 
backbone does not lose its force indeterminacy at low 
pressure, even in the slow compression limit, except for 
fi — > +00. The degree of force indeterminacy, H, de- 
creases to about 14% of the number of backbone degrees 



of freedom in poorly coordinated systems with fi = 0.3 
On computing H it is necessary to take into account 
the contribution of divalent grains, which define local- 
ized (harmless) floppy modes. Consequently the value of 
the backbone coordination number z* corresponding to 
H = is slightly smaller than 4. 



V. DISCUSSION AND PERSPECTIVES 

The most important novel feature of our simulation re- 
sults is the wide variability of coordination numbers for 
the same density in isotropic packings. Most often, dense 
numerical samples are obtained on suppressing friction: 
A-type packings are simulated, in which the coordination 
number is high. C-type systems have about the same 
density, but a coordination number as low as in the loos- 
est states D obtained by direct compression. In order 
to simulate dense laboratory samples, should one use A 
or C conflgurations ? The answer depends of course on 
how laboratory samples are assembled. Our results show 
that vibrated ones are likely to have smaller coordina- 
tion numbers than lubricated ones for the same density. 
In paper III, we compare elastic properties of our B and C 
states to the ones measured by Jia and Mills [sl] on glass 
bead samples either compacted by shaking the container 
or assembled with a lubricant. 

In the X-ray tomography experiments by Aste et 
al. 0], sphere packings are imaged with a resolution 
(voxel size) of about 4% of nominal diameter a, while 
the diameter distribution extends at least to ±0.03a. In 
spite of serious efforts to eliminate the influence of size 
distribution by deconvoluting correlation data, their esti- 
mation of coordination numbers are well above the upper 
limit 6 in dense samples, which is in principle impossible 
under a low pressure. It seems that such experiments 
only provide access to the largest of the three length 
scales mentioned in Section IIV E 21 and are thus unable 
to distinguish between A-like and C-like microstructures. 
We shall see in paper III that measurements of elas- 
tic moduli are much better suited to obtain information 
on coordination numbers by experimental means. 

However, it is flrst necessary to assess the influence of 
the pressure level on the packing inner states. As hinted 
by the results of Table IIIIl a quasistatic compression af- 
fects the force distribution and the level of friction mobi- 
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lization. Elastic properties being usually measured above 
a certain confinement level (typically, a few tens of kPa) , 
the necessary study of the effects of a quasistatic com- 
pression is carried out in paper II [37l|. 

Beyond the elastic properties, which characterize the 
response to small load increments, the quasistatic, elasto- 
plastic mechanical behavior of packings prepared with 
different microstructures should also be studied. Peak 
deviator strength and dilatancy normally correlate with 
initial packing density dH, [S [11 [M], [23] . But for one 
given density, what is the influence of the coordination 
number on the stress-strain curves ? 

Granular systems are often packed under gravity by 
pouring samples in containers, and such processes, which 
do not necessarily produce homogeneous states [s^ 
should be studied by numerical simulations too, and 
the analysis of dynamical effects in the assembling stage 
should be pursued. Adhesive contact forces, as in wet 
granular assemblies ^] , can also greatly affect the prepa- 
ration of sohd granular packings [99|. Other perspec- 
tives to the present work are the investigation of the mi- 
crostructure of polydispe rse s ystems, and of assemblies 
of non-spherical particles jlOOl | . 



APPENDIX A: CONTACT ELASTICITY AND 
FRICTION 

The contact law between spherical elastic bodies, with 
a Coulomb criterion for friction applied locally (to the 
surface force densities), leads to complicated history- 
dependent force-displacement relationships [s^, Il0l| . 
Even in some cases with no slip anywhere in the con- 
tact region, the tangential stiffness Kt of a contact was 
shown [102| | to depend on the past history of the contact 
loading, and to change according to the direction of the 
displacement increment. Strictly speaking, the response 
of intergranular contacts, even to arbitrary small load 
increments, should not be termed "elastic" . The sim- 
plified law we adopted involves a tangential stiffness Kt 
depending on the normal deflection h, but independent 
of the current mobilization of friction. This is the same 
approximation as used in [2^, [s^l : the value of Kt is the 
correct one in the absence of elastic relative tangential 
displacement, when T = 0. 

However, as stressed in [i^l, such a model is thermo- 
dynamically inconsistent, for the elastic energy might in- 
crease at no cost. Consider, e.g., quasistatically reduc- 
ing h at constant Ju^, thereby, according to this contact 
model, reducing normal force N at constant T, without 
reaching the Coulomb limit. The recoverable elastic en- 
ergy stored in the contact is given by 



2 1 

5 2 Kt 



(Al) 



which grows as Kt decreases, without the external force 
doing any work, thus implying a net creation of energy. 
To avoid such effects, T is rescaled (as advocated by O. 



Walton [1Q3[), whenever N decreases to iV — AiV, down 
Kt{N - l\N) 

to T ; — , before accounting for tangential rela- 

Kt{N) ^ ^ 

five displacement increments. No such rule applies to in- 
creasing normal force cases. Such a procedure was shown 
by Elata and Berryman 40] to systematically produce 
energy dissipation in cyclic loadings of the contact. 

Such peculiarities of the contact law affect the form 
and, in fact, the very definition of an elastic response of 
the contact network, an issue which will be discussed in 
paper III. 



APPENDIX B: TRANSPORT OF CONTACT 
FORCES DUE TO PARTICLE MOTION 

In molecular dynamics calculations, as well as in static 
approaches (as outlined in Section lll G\ one has to relate 
small contact force increments in any contact to the small 
displacements and rotations A^^ of the grains. 

Increments A(iVijnij) and AT^j of the normal and tan- 
gential parts of the force in the contact between grains 
i and j have two different origins: they stem from the 
contact law, as written down in Section |TT1 and also from 
the motion of the particle pair. As the grains move, so 
docs the deformed contact region, and therefore the re- 
sulting contact force changes. The relevant formulae are 
derived and written below for small increments in the 
static case. For dynamical computations, displacements 
are to be replaced by velocities, and increments by time 
derivatives. 

The normal force variation is simply 



A(iVyny) = AiVy-n,, 



iVy-Ariy, 



(Bl) 



with 



An, 



(i - ® n.y) • (uj - u — I • Tij) , (B2) 



while AA^ij is related by the Hertz law to the variation 
in the normal deflection of the contact. 

For the tangential reaction we introduce the decompo- 
sition 



AT,, 



AT 



(1) 



AT 



(2) 



Increments with superscript (1) are associated, via the 
contact law, to the relative displacement of the contact 
point, which defines the constitutive part of the stiffness 
matrix discussed in Section III CI and superscript (2) la- 
bels increments of kinematic origin. We assume that the 
magnitude of the contact force is unchanged in the ab- 
sence of relative displacement at the contact {Suij — 0), 
and thus we write 



AT 



.(2) 



: A% X Ty, 



(B3) 



A0y denoting the rotation of the contact region. A6'ij 
can be split in a rolling part A^^J^-*, orthogonal to n^ . 
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and a pivoting one, A^^-^ , along riy. A9lj ' is deter- 
mined by the incremental change in 11^^: 

M\f = X An,,. (B4) 

As to the pivoting part A6'j^^ , it is natural to equate it 
to the average rotation of the two particles around the 
normal direction: 

A0(f = in,, .(A^. + A^,)!!,,-. (B5) 

This choice is such that the rotation of the contact force 
coincides with the rotation of the pair in contact if both 
objects move together like one rigid body (a condition of 
objectivity [49]). 

Injecting (|B^ into Eqns. (|B4)) and (|B5l) . one readily 
obtains the appropriate formula for AT^y , 

AT^f =-[T,,.(u,-u,-e.r,,)]^ 

+ i[(A0, + A^^,)-n,,](ny x T„). 

With the notations of Section III C[ the contribution of 
contact i,j to the 3x3 diagonal block of K*-^-* which 
expresses the dependence of F"''* on displacement is 
the non-symmetric tensor 

In general, the geometric stiffness matrix K'^' is thus not 
symmetric, except in frictionless sphere packings, which 
are analogous to central force networks. We also note 
that terms of order N/R or ||T||/i? in K'^' correspond 

to terms of order Kn or Kt in K'^'', which are always 
very much larger. 

In the frictionless case, K^^-* is a symmetric, negative 
matrix if forces are repulsive, as discussed by Alexan- 
der [9l|. Any mechanism on the backbone leads to an 
instability: the potential energy of the externally applied 
load is strictly decreasing in that motion. This destabiliz- 
ing effect can also be directly established in the rigid case, 
as shown in [3^ . This is the reason why stable packings 
of frictionless spheres in equilibrium under some exter- 
nally applied load are devoid of mechanisms involving 
the backbone. 

In general stiffness matrices were discussed by Kuhn 
and Chang [i^, and by Bagi [s^. Those authors gave 
general results for K*-^-* with particles of arbitrary shapes, 
which coincide with ours in the case of spherical balls. 



APPENDIX C: STRESS-CONTROLLED 
MOLECULAR DYNAMICS 

It might be noted that the original Parrinello- Rahman 
equations slightly differ from ours. First, Eqn. ^ is writ- 



ten down with an additional term 

m,sf) = -^^;(")-2^^m,i("). (CI) 

Then, eqn. ([7]) is written with a different stress tensor, 
TT, the definition of which involves a particular reference 
value Lq of the cell dimensions, for which the volume is 
Hq. tt is related to the Cauchy stress tensor S by 

TT is a symmetric tens or k nown as the second Piola- 
Kirchhoff stress tensor Il0 4| (al so ca lled thermodynamic 
tension by some authors [lOSl Il06j |). Tensor tt can be 
used to express the power V /VLq of internal forces, per 
unit volume in the reference (undeformcd) configuration. 

The metric tensor Q_ — '^Lq ^ • "^L • L • Lg^^ expresses 
the distance between current points as a function of their 
coordinates in the reference configuration. The differ- 
ence between Q and the unit tensor defines the Green- 
Lagrange strain tensor |l04l |. e, as 

e = -^(£-l). (C3) 

Then tt is such that, for whatever strain history 

V^nog: |. (C4) 
fL")^ 

If the last term of Eqn. ^ is replaced by -—2 — , and 

if (jCip is used instead of ©, then, in the case when 
interparticle forces derive from a potential V, function 
of particle positions and orientations, the system of 
equations conserves the total energy 

(C5) 

Such equations would tend to impose a constant Piola- 
Kirchhoff stress. 

Granular assemblies are however dissipative systems, 
and energy conservation is not a crucial issue as in molec- 
ular systems. In practice, we observed that omission of 
the extra term of (jCip as well as control of Cauchy, rather 
than Piola-Kirchhoff stresses, did not affect the approach 
to equilibrium configurations. Yet, on considering small 
motions and elastic properties close to an equilibrium 
configuration, one may prefer dealing with external forces 
deriving from a potential. We note that this is indeed the 
case if we use Eqns. (|Cip and ([7]) with an isotropic S, 
S = PI, in which case the external stress control is asso- 
ciated with potential energy PVl (instead of the last term 
in Eqn.[C5|. 
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APPENDIX D: COMPARISON OF LOCAL BOND 
ORDER PARAMETERS WITH EXPERIMENTAL 
OBSERVATIONS 

Although the samples were not isotropic, owing to the 
role of gravity in the preparation stage, and not perfectly 
homogeneous, because of lateral walls, the X-ray tomog- 
raphy experiments by Aste et al. [7[ provided an un- 
precedented wealth of results on the geometry of sphere 
packings. The local arrangement of neighbors around 
one bead were classified according to the values taken 
by the pair {Q^^'^ {i) , Q^^^ (i) , for different choices of the 
distance dc that defines the bonds. For each sample 
and choice of dc, the proportion of spheres having the 
most frequent range of values {Q4 ± 0.05, Qe ± 0.05) cor- 
responding to some typical disordered arrangement was 
recorded, as well as the frequency of occurrence of values 
(0.191±0.05, 0.574±0.05) and (0.097±0.05, 0.485±0.05) 
respectively corresponding to fcc-like and hcp-like local 
ordering. Those values are compared to the ones we 
observed in numerical samples of similar density in ta- 
bles [Vlll and |VIIll (we kept the same {QA^Qe) couple). 

TABLE VIL Most frequently values {Q4{i), Qeii)) observed in 
the measurements of [J in a dense sample ($ — 0.640 ±0.005, 
called sample F in 0]), and fraction (dis.) of local configura- 
tions within interval (Q4 ±0.05, Qe ±0.05) obtained in exper- 
iments and in numerical simulations with similar densities : 



A, A' and C. 



Sample 


dc 


(Q4, Qe,) 


dis (%) 


fee (%) 


hep (%) 


Aste et al. 


1.1 


(0.23, 0.44) 


31 


6 


4 


dense 


1.2 


(016, 0.45) 


38 


4 


12 


$ ~ 0.640 


1.3 


(0.13, 0.42) 


43 


1 


17 


(F in [7]) 


1.4 


(0.10, 0.38) 


47 


3 


13 


A 


1.1 


(same 


32 


< 10"^ 


4 


A 


1.2 


values 


37 


< 10-^ 


5 


A 


1.3 


as 


43 


< 10"^ 


8 


A 


1.4 


above) 


46 


< 10^2 


12 


A' 


1.1 


(same 


31 


< 10"^ 


5 


A' 


1.2 


values 


40 


< 10"^ 


6 


A' 


1.3 


as 


45 


< 10-2 


10 


A' 


1.4 


above) 


48 


< 10-2 


15 


C 


1.1 


(same 


31 


< 10-2 


4 


C 


1.2 


values 


37 


< 10-2 


5 


c 


1.3 


as 


43 


< 10-2 


8 


c 


1.4 


above) 


46 


< 10-2 


13 



The fraction of beads with the typical disordered config- 
uration of neighbors (marked dis. in the tables) are very 
close in numerical packings and in the experimental one 
of the same densities, and the frequency of occurrence of 
hcp-like local environments also compares well, although 
it does not seem to share the same dependence on the 
threshold distance dc defining bonds. However, the small 
fraction of fcc-like beads observed in the laboratory is 



TABLE VIIL Same as table IVnl for experimental samples of 
lower densities, one with "I> — 0.626 ± 0.008 (called sample 
D in [3]), another with $ = 0.596 ± 0.006 (called sample B 
in 01), to which values obtained in simulated samples B and 
D, of similar densities, are respectively compared. 



Sample 


dc 


{Q4, Qg) 


dis (%) 


fee (%) 


hep (%) 


Aste et al. 


1.1 


(0.25, 0.44) 


28 


4 


1 




1.2 


(0.19, 0.44) 


35 


2 


7 


<E> ~ 0.626 


1.3 


(0.15, 0.40) 


42 


1 


11 


(D in [7]) 


1.4 


(0.11, 0.36) 


46 


1 


8 


B 


1.1 


(same 


30 


< 10-2 


2.7 


B 


1.2 


values 


38 


< 10-2 


7.6 


B 


1.3 


as 


43 


< 10-2 


10 


B 


1.4 


above) 


46 


< 10-2 


10 


Aste et al. 


1.1 


(0.30, 0.45) 


24 


3 


1 


loose 


1.2 


(0.23, 0.44) 


32 


2 


3 


$ ~ 0.596 


1.3 


(0.16, 0.38) 


37 


1 


5 


(B in [7]) 


1.4 


(0.14, 0.35) 


43 


2 


5 


D 


1.1 


(same 


24 


< 10-2 


1.0 


D 


1.2 


values 


32 


< 10-2 


4.4 


D 


1.3 


as 


37 


< 10-2 


6.1 


D 


1.4 


above) 


43 


< 10-2 


7.5 



absent in the simulations. Many cirmcumstances can be 
invoked to explain these differences, including of course 
the different packing history of the numerical and ex- 
perimental samples, which in the the latter case involves 
gravity and anisotropy. It can be remarked once again 
that A' samples are a little more ordered than A ones 
(with slightly larger fractions of hcp-like local configu- 
rations), from which C samples are quite indistinguish- 
able, as the quantities measured here do not depend on 
whether pairs of neighbors are actually in contact. 



APPENDIX E: ANALYSIS OF THE FREE 
MOTION OF DIVALENT SPHERES 

We first give here the appropriate formulae to describe 
the free motion of divalent grains, then report on numer- 
ical stability tests. 

The equations are specialized to the case of equal-sized 
spheres of dimater a, as in all the simulations reported in 
the present paper. Let i denote the label of the divalent 
grain in contact with its neighbors labelled j and k. The 
line joining the centers of j and k is parallel to unit vector 
e, defined as 

\\nik — UijW ' 

and the distance D of the center of i to this line is 
D = ay/l - (e • Uij^. 
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ujo denoting the angular velocity of the center of i about 
this axis, the translational velocity of i, in its free motion 
depicted on Fig. [TOl will be 

V, = cJoDt, (El) 

the unit vector t being orthogonal to the plane containing 
the centers of i, j, and k: 

^ ^ n,,j X n,fc 
\\nik X riyll ' 

Its angular velocity will be 

= 2LUoe. (E2) 

With such a choice, the instantaneous velocity of the (ma- 
terial) contact points between i and j, or between i and 
k satisfy: 

v,j + S2i X -Uij = v,j + S2j X -rijfc = 0, 
2 2 

as requested in a relative motion which is a combination 
of rolling and pivoting. 

It is easy to check that the rules defined in Appendix [Bl 
ensure that the tangential components of the contacts i, j 
and i, k rotate with the contact points around the axis 
joining the centers of k and j with angular velocity luq. 
In other words, the geometric stiffness matrix K*-^-* does 
not determine whether the mechanism associated with a 
two-coordinated bead is stable. 

We check for stability with numerical means, as fol- 
lows. Starting from an equilibrium configuration, we first 
choose the potentially dangerous mechanisms, those in- 
volving relatively large contact forces, of the order of the 
average normal force or even larger. Thus, grains j and 
k undergo significant changes in the moment of the con- 
tact force with the mobile grain i. Then one such mobile 
divalent bead is attributed a velocity and a angular ve- 
locity according to Egns. [ET] and lE2l with a value of ojo 
small enough for the centrifugal acceleration to be negli- 
gible (equal to 10~*||Fij ||). The evolution of the whole 
packing under constant stress is then simulated with 
MD. Such numerical experiments were performed with 
the most fragile packings, D samples under low confining 
pressures. In all cases studied, as expected, the motion 
of the mechanism entails very slow changes in the con- 
figuration, if any. The mobile grain maintains a constant 
angular velocity while nearly exactly following its circu- 
lar trajectory for a long time, with hardly any change in 
kinetic energy. These calculations are rather slow and 
costly, and we therefore limited our investigations to 10 
tests for 2 pressure levels in the D series, P = 1 kPa 
and P = 10 kPa. At the lowest pressure P — 1 kPa, 



corresponding to k ~ 181000, these motions were often 
observed to lead to a small rearrangement of the packing, 
in which kinetic energy spreads over all degrees of free- 
dom, a significant fraction of the contacts, up to 25%, go 
through a sliding stage, the contact network is slightly 
modified and the system restabilizes in a slightly differ- 
ent configuration, with a small density increase (typically 
of order 10~^). In other cases, the freely moving grain 
stops when it collides with a third grain other than the 
two with which it is maintaining contact. The system 
then finds a new equilibrium configuration without rear- 
ranging, only a few contacts temporarily reach a sliding 
status. On repeating similar tests at a larger confining 
pressure, P = 10 kPa (still close to the rigid limit), the 
occurrence of this second scenario became much more 
frequent than the first, which was never observed in the 
10 tests performed. 

The difference between stable and unstable cases is 
better appreciated on redoing the tests with a modified 
MD calculation method, in which only the initially ex- 
isting contacts are taken into account. Thus one only 
investigates the properties of the pre-existing contact net- 
work. If it breaks, the system globally falls apart, and 
nearly all contacts in the packing eventually open. This 
is the unstable case. In the stable case the mobile par- 
ticle can turn several times around the line of centers 
of its two contacting neighbors without notable changes 
in kinetic energy and the contact network is maintained. 
This behavior is illustrated on Fig. [121 which displays 
trajectories of mobile grains in the plane orthogonal to 
e. Note that such a procedure reveal instabilities that 
are prevented by the appearance of a third contact of the 
mobile grain, and hence overestimates the frequency of 
occurrence of unstable configurations. 8 out of 10 such 
tests led to an instability in D samples under 1 kPa. This 
proportion fell to 2 out of 10 under 10 kPa. 

Two possible conclusions may be drawn. On the one 
hand we may deem the D configurations under low pres- 
sures imperfectly stabilized, as some free motions might 
eventually cause configurational changes. Or, since any- 
way the evolution is so slow, it may be pointed out, on 
the other hand, that the slightest amount of dissipation 
in rolling would stop the free motion. In realistic systems 
the velocity a a free rolling motion always decays in time. 

Since the obtention of stable, equilibrated states 
in which all divalent grains have been made three- 
coordinated would be computationnally very costly, and 
as the occurrence of small instabilities related to such 
mechanisms appear to decrease fastly under growing con- 
finement, we adopted the second attitude and regarded 
D and C configurations with a few percent of divalent 
particles as acceptable equilibrium states. 
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FIG. 12: (Color online) Examples of trajectories of mobile 
divalent grains in the plane orthogonal to vector e. All tra- 
jectories begin a,t x — D, y = 0. Continuous lines depict 
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up, in 2 cases out of the 4 represented, in instabilities (black 
lines), while complete circles are observed in the other, stable 
ones (red). Dots mark the same trajectories, which are in 
practice arrested by other contacts when contact creation is 
taken into account. 



Granular Media (Wiley- VCH, Berlin, 2004). 
[3] R. Garcia Rojo, H. J. Herrmann, and S. McNamara, 
eds.. Powders and Grains 2005 (Balkema, Leiden, 
2005). 

[4] A. J. Liu and S. R. Nagel, eds.. Jamming and rheology 

(Taylor & Francis, New York, 2001). 
[5] P. Richard, P. Philippe, F. Barbe, S. Bourles, 

X. Thibault, and D. Bideau, PRE 68, 020301(R) (2003). 
[6] T. Aste, M. Saadatfar, A. Sakellariou, and T. J. Senden, 

Physica A 339, 16 (2004). 
[7] T. Aste, M. Saadatfar, and T. J. Senden, PRE 71, 

061302 (2005). 

[8] D. Cumberland and R. Crawford, The Packing of Par- 
ticles (Elsevier, Amsterdam, 1987). 
[9] D. Bideau and J. Dodds, eds.. Physics of Granular Me- 
dia (Nova Science Publishers, 1991). 

[10] S. Torquato, T. M. Truskett, and P. G. Debenedetti, 
Physical Review Letters 84, 2064 (2000). 

[11] D. M. Wood, Soil Behaviour and Critical State Soil Me- 
chanics (Cambridge University Press, 1990). 

[12] J. Biarez and P.-Y. Hicher, Elementary Mechanics of 
Soil Behaviour (A. A. Balkema, Rotterdam, 1993). 

[13] J. K. Mitchell, Fundamentals of soil behavior (Wiley, 
New York, 1993). 

[14] M. Oda, Soils and Foundations 12, 17 (1972). 

[15] J. R. F. Arthur and B. K. Menzies, Geotechnique 22, 
115 (1972). 

[16] S. Fukushima and F. Tatsuoka, Soils and Foundations 



24, 40 (1984). 

[17] L. Vanel, D. Howell, D. Clark, R. P. Behringer, and 
E. Clement, Phys. Rev. E 60, R5040 (1999). 

[18] P. A. Cundall and O. D. L. Strack, Geotechnique 29, 
47 (1979). 

[19] R. J. Bathurst and L. Rothenburg, Mechanics of Mate- 
rials 9, 65 (1990). 

[20] F. Radjai and S. Roux, in [ipl, pp. 21-24. 

[21] F. Radjai and S. Roux, in [37pP- 165-187. 

[22] L. Rothenburg and N. P. Kruyt, International Journal 
of Solids and Structures 41, 5763 (2004). 

[23] J.-N. Roux and F. Chevoir, Bulletin des Laboratoires 
des Fonts et Chaussees 254, 109 (2005). 

[24] L. E. Silbert, D. Erta§, G. S. Grest, T. C. Halsey, and 
D. Levine, Physical Review E 65, 031304 (2002). 

[25] L. E. Silbert, G. S. Grest, and J. W. Landry, Physical 
Review E 66, 061303 (2002). 

[26] H. Makse, D. Johnson, and L. Schwartz, Physical Re- 
view Letters 84, 4160 (2000). 

[27] C. Thornton, Geotechnique 50, 43 (2000). 

[28] A. S. J. Suiker and N. A. Fleck, ASME Journal of Ap- 
plied Mechanics 71, 350 (2004). 

[29] H. Makse, N. Gland, D. Johnson, and L. Schwartz, 
Physical Review Letters 83, 5070 (1999). 

[30] H. A. Makse, N. Gland, D. L. Johnson, and L. Schwartz, 
Physical Review E 70, 061302 (2004). 

[31] C. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, 
Physical Review E 68, 011306 (2003). 

[32] J.-N. Roux, Physical Review E 61, 6802 (2000). 

[33] N. S. Rad and M. T. Tumay, ASTM Journal of Geotech- 
nical Testing 10, 31 (1987). 

[34] N. Benahmed, J. Canou, and J.-C. Dupla, Comptes- 
Rendus Academie des Sciences, Mecanique 332, 887 
(2004). 

[35] E. R. Nowak, J. B. Knight, E. Ben-Naim, H. M. Jaeger, 
and S. R. Nagel, Physical Review E 57, 1971 (1998). 

[36] P. Philippe and D. Bideau, Europhysics Letters 60, 677 
(2002). 

[37] I. Agnolin and J.-N. Roux, Internal states of model 
isotropic granular packings: compression and pressure 
cycles, first companion paper, second in the series 
(2007). 

[38] I. Agnolin and J.-N. Roux, Internal states of model 
isotropic granular packings: elastic properties, second 
companion paper, third in the series (2007). 

[39] K. L. Johnson, Contact Mechanics (Cambridge Univer- 
sity Press, 1985). 

[40] D. Elata and J. G. Berryman, Mechanics of Materials 
24, 229 (1996). 

[41] E. Somfai, J.-N. Roux, J. Snoeijer, M. van Hecke, and 
W. van Saarloos, PRE 72, 021301 (2005). 

[42] M. Parrinello and A. Rahman, Journal of Applied 
Physics 52, 7182 (1981). 

[43] J. Christoffersen, M. M. Mehrabadi, and S. Nemat- 
Nasser, Journal of Applied Mechanics 48, 339 (1981). 

[44] J. -P. Hansen and I. R. McDonald, Theory of Simple Liq- 
uids (Academic Press, 1986). 

[45] P. Cundall, in Micromechanics of Granular Materials, 
edited by M. Satake and J. T. Jenkins (Elsevier, Ams- 
terdam, 1988), pp. 113-123. 

[46] B. Servatius and H. Servatius, in [sil], pp. 1-19. 

[47] S. McNamara, R. Garcia Rojo, and H. J. Herrmann, 
Physical Review E 72, 021304 (2005). 

[48] S. McNamara and H. J. Herrmann, Physical Review E 



29 



74, 061303 (2006). 
[49] M. R. Kuhn and C. S. Chang, International Journal of [77 
Solids and Structures 43, 6026 (2006). [78 
K. Bagi, Granular Matter 9, 109 (2007). [79; 
M. F. Thorpe and P. M. Duxbury, eds.. Rigidity The- 
ory and Applications, Fundamental Materials Research [80 
(Kluwer Academic, 1998). 

A. Donev, S. Torquato, and F. H. Stillinger, PRE 71, 
011105 (2005). [81 
G. Combe and J.-N. Roux, in Deformation charac- 
teristics of geomaterials, edited by H. di Benedetto, [82 
T. Doanh, H. Geoffroy, and C. Sauzeat (Swets and 
Zeithnger, Lisse, 2003), pp. 1071-1078. [83; 

F. Tatsuoka, M. Sakamoto, T. Kawamura, and 
S. Fukushima, Soils and Foundations 26, 65 (1986). [84; 

G. Reydellet and E. Clement, Physical Review Letters [85 
86, 3308 (2001). 

G. Combe, Mecanique des materiaux graulaires et orig- [86 
ines micro scopiques de la deformation (Presses du Lab- 
oratoire Central des Fonts et Chaussees, Paris, 2002). [87 
GdR MIDI, European Physical Journal E 14, 341 
(2004). [ 
F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and 
F. Chevoir, Physical Review E 72, 021309 (2005). [89; 

0. Pouliquen, C. Cassar, Y. Forterre, P. Jop, and 
M. Nicolas, in Q, pp. 859-865. 

P. Jop, Y. Forterre, and O. Pouliquen, Journal of Fluid [90 
Mechanics 541, 167 (2005). 

P. Jop, Y. Forterre, and O. Pouhquen, Nature 44, 727 [91 
(2006). [92' 

B. J. Alder and T. E. Wainwright, J. Chem. Phys. 31, [93' 
459 (1959). 

B. D. Lubachevsky and F. H. Stillinger, Journal of Sta- [94 
tistical Physics 60, 561 (1990). [95; 

B. D. Lubachevsky, F. H. Stillinger, and E. N. Pinson, 
Journal of Statistical Physics 64, 501 (1991). [96; 
A. Donev, S. Torquato, F. H. Stillinger, and R. Con- 
nelly, Journal of Computational Physics 197, 139 
(2004). [97 

1. Volkov, M. Cieplak, J. Koplik, and J. R. Banavar, 
Physical Review E 66, 061401 (2002). [l 
S. Emam, J.-N. Roux, J. Canou, A. Corfdir, and J.-C. 
Dupla, in 0], pp. 49-52. [99; 
N. Benahmed, Ph.D. thesis, Ecole Nationale des Fonts 
et Chaussees, Champs-sur-Marne, France (2001). [100 
V. Luchnikov, A. Gervois, P. Richard, L. Oger, and J.- 
P. Troadec, Journal of Molecular Liquids 96-97, 185 [101 
(2002). 

R. Connelly, Structural Topology 14, 43 (1988). [102 
A. R. Kansal, S. Torquato, and F. H. Stillinger, Physical 
Review E 66, 041109 (2002). [103; 
S. Ouaguenouni and J.-N. Roux, Europhysics Letters, [104 
39, 117 (1997). 

J.-N. Roux, in Proceedings of the Saint-Venant Sym- [105 
posium on Multiple Scale Analysis and Coupled Physi- [106 
cal Systems (Presses de I'Ecole Nationale des Fonts et 
Chaussees, Paris, 1997), pp. 577-584. [107 

C. F. Moukarzel, Physical Review Letters 81, 1634 
(1998). [108; 
C. F. Moukarzel, in [H, pp. 125-142. 
A. Tkachenko and T. A. Witten, Physical Review E 60, 



627 (1999). 

C. F. Moukarzel, Granular Matter 3, 41 (2001). 

C. F. Moukarzel, in f^, pp. 23-43. 

L. E. Silbert, A. J. Liu, and S. R. Nagel, Physical Re- 
view E 73, 041304 (2006). 

Determination des masses volumiques minimales et 
maximales des sols non coherents, French Standard 
NFP 94-059 (1992). 

G. Y. Onoda and E. G. Liniger, Physical Review Letters 
64, 2727 (1990). 

M. M. Kohonen, D. Geromichalos, M. Scheel, C. Schier, 
and S. Herminghaus, Physica A 39, 7 (2004). 

H. P. Zhang and H. A. Makse, Physical Review E 72, 
011301 (2005). 

X. Jia and P. Mills, in [lo3], pp. 105-112. 

G. T. Nolan and P. E. Kavanagh, Powder Technology 

72, 149 (1992). 

W. M. Visscher and M. Bolsterli, Nature 239, 504 
(1972). 

R. Jullien and P. Meakin, Europhysics Letters 4(12), 
1385 (1987). 

G. C. Barker and A. Mehta, Physical Review A 45, 3435 
(1992). 

A. Mehta, G. C. Barker, and J. M. Luck, Journal of Sta- 
tistical Mechanics: Theory and Experiment p. P10014 
(2004). 

I. Bratberg, F. Radjai, and A. Hansen, Physical Review 
E 66, 031303 (2002). 

S. Alexander, Phys. Rep. 296, 65 (1998). 

S. Ouaguenouni and J.-N. Roux, in [10^, pp. 188-191. 

F. Radjai, M. Jean, J. -J. Moreau, and S. Roux, Phys. 

Rev. Lett. 27, 274 (1996). 

S. J. Antony, Phys. Rev. E 63, 011302 (2000). 

D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Phys. Rev. 
E 57, 3164 (1998). 

D. L. Blair, N. W. Mueggenburg, A. H. Marshall, 

H. Jaeger, and S. R. Nagel, Physical Review E 63, 
041304 (2001). 

F. Radj ai, D. Wolf, S. Roux, M. Jean, and J.-J. Moreau, 
in [loi, pp. 169-179. 

F. Radjai, D. E. Wolf, M. Jean, and J.-J. Moreau, Phys- 
ical Review Letters 80, 61 (1998). 

F. A. Gilabert, J.-N. Roux, and A. Castellanos, Phys. 
Rev. E 75, 011303 (2007). 

A. Donev, R. Connelly, F. H. Stillinger, and 
S. Torquato, PRE 75, 051304 (2007). 
R. D. Mindlin and H. Deresiewicz, ASME Journal of 
Applied Mechanics 20, 327 (1953). 

D. L. Johnson and A. N. Norris, Journal of Mechanics 

and Physics of Solids 45, 1025 (1997). 

O. Walton, private communication. 

J. Salengon, Handbook of continuum mechanics. General 

concepts, Thermoelasticity. (Springer, 2001). 

J. R. Ray and A. Rahman, JCP 76, 4423 (1984). 

A. Lemaitre and C. Maloney, Journal of Statistical 

Physics 123, 415 (2006). 

Y. Kishino, ed.. Powders and Grains 2001 (Swets & 
Zeitlinger, Lisse, 2001). 

D. Wolf and P. Grassberger, eds.. Friction, Arching, 
Contact Dynamics (World Scientific, Singapore, 1997). 



